Autor: Szymon Pawłowski

Naszym zadaniem jest przeanalizowanie zbioru danych medycznych oraz stworzenia modelu klasyfikującego potencjalnego pacjenta ze względu na fakt wystąpienia udaru mózgu.

Udar mózgu to zespół objawów klinicznych związanych z nagłym wystąpieniem ogniskowego lub uogólnionego zaburzenia czynności mózgu, powstały w wyniku zaburzenia krążenia mózgowego i utrzymującego się ponad 24 godziny.

Spis treści

• Biblioteki i funkcje
• Dane
• EDA
  • stroke
  • id
  • gender
  • age
  • hypertension
  • heart_disease
  • ever_married
  • work_type
  • residence_type
  • avg_glucose_level
  • bmi
  • smoking_status
• Inżynieria cech
  • bmi
  • gender
  • age
  • ever_married
  • work_type
  • residence_type
  • avg_glucose_level
  • smoking_status
• Analiza korelacji
• Polynomial Features
• Balansowanie zbioru
  • SMOTE
  • SMOTE+Tomek
  • SMOTE+ENN
• Metryki
• Modelowanie
  • Wybór techniki balansowania
  • Selekcja zmiennych przy użyciu LASSO
  • Hyperparameters tuning
  • Wybór thresholdu $p$
• Przykładowy pipeline

Biblioteki i funkcje

powrót do spisu treści

import pandas as pd
import numpy as np
import joblib

from plotly.subplots import make_subplots
import plotly.graph_objects as go
import plotly.express as px
from matplotlib import pyplot as plt
import seaborn as sns

from sklearn import preprocessing
from sklearn.utils import shuffle
from sklearn.model_selection import train_test_split
from imblearn.over_sampling import SMOTE
from imblearn.over_sampling import SMOTENC
from imblearn.over_sampling import ADASYN
from imblearn.combine import SMOTETomek
from imblearn.combine import SMOTEENN
from imblearn.under_sampling import EditedNearestNeighbours
from sklearn.pipeline import Pipeline
from sklearn import metrics
from sklearn.model_selection import GridSearchCV

from sklearn.linear_model import LogisticRegression
from sklearn.svm import SVC
from sklearn.feature_selection import SequentialFeatureSelector
from sklearn.ensemble import RandomForestClassifier
from sklearn.feature_selection import RFE
from sklearn.model_selection import cross_validate
from sklearn.feature_selection import SelectFromModel

from yellowbrick import classifier
from yellowbrick.classifier import ClassificationReport
from yellowbrick.classifier.rocauc import roc_auc

import warnings
warnings.filterwarnings("ignore")

def plotting(data: pd.DataFrame, x: str, type: str, y=None, cl = 'stroke'):
  """
    Function for plotting data in specific way

  Args:
    data - dataset to plot

    x - column for x-axis
    
    y - column for y-axis
    
    type - type of plot 
      values: {'hist', 'box', 'hist_box', 'class_hist', 'class_box', 'class_hist_box', 'scatter', 'heatmap', 'lineplot'} 

    cl - class

  Returns:
    Plot in plotly
  """

  def plot_class_hist():
    """
      Plots histogram with stroke classes in plotly
    """
    
    fig = go.Figure()
    fig.add_trace(go.Histogram(x = data.loc[data[cl]==0, x], name='No '+cl))
    fig.add_trace(go.Histogram(x = data.loc[data[cl]==1, x], name=cl))

    fig.update_layout(barmode='overlay')
    fig.update_traces(opacity=0.75)
    
    return fig
  
  def plot_class_box():
    """
      Plots boxplot with stroke classes in plotly
    """
    
    fig = go.Figure()
    fig.add_trace(go.Box(y = data.loc[data[cl]==0, x], name='No '+cl))
    fig.add_trace(go.Box(y = data.loc[data[cl]==1, x], name=cl))

    fig.update_layout(barmode='overlay')
    fig.update_traces(opacity=0.75)
    
    return fig
    

  if type=='hist':
    fig = px.histogram(data, x=x)
    fig.show()

  elif type=='box':
    fig = px.box(data, x=x)
    fig.show()
  
  elif type=='scatter':
    data['temp'] = data[cl].astype('str')

    fig = px.scatter(data, x=x, y=y, color='temp')
    fig.update_layout(title_text="Scatterplot for " + x)
    fig.update_traces(marker_size=10)
    
    fig.show()
  
  elif type=='heatmap':
    corr = data.corr()
    f, ax = plt.subplots(figsize=(16, 12))
    mask = np.triu(np.ones_like(corr, dtype=bool))
    cmap = sns.diverging_palette(230, 20, as_cmap=True)
    sns.heatmap(corr, annot=True, mask = mask, cmap=cmap)

  elif type=='hist_box':
    fig = make_subplots(rows=1, cols=2)

    fig.add_trace(go.Histogram(x = data[x], name=x), row=1,col=1)
    fig.add_trace(go.Box(y = data[x], name=x), row=1, col=2)
    
    fig.update_layout(barmode="overlay")
    fig.update_traces(opacity=0.7, row=1, col=1)
    fig.update_layout(title_text="Histogram and Boxplot for " + x)
    fig.show()

  elif type=='class_hist':
    
    fig = plot_class_hist()
    fig.show()
  
  elif type=='class_box':
    
    fig = plot_class_box()
    fig.show()

  elif type=='class_hist_box':
    fig = make_subplots(rows=1, cols=2)

    fig.add_trace(go.Histogram(x = data.loc[data[cl]==0, x], name='No '+cl), row=1,col=1)
    fig.add_trace(go.Histogram(x = data.loc[data[cl]==1, x], name=cl), row=1, col=1)
    fig.add_trace(go.Box(y = data.loc[data[cl]==0, x], name='No '+cl), row=1, col=2)
    fig.add_trace(go.Box(y = data.loc[data[cl]==1, x], name=cl), row=1, col=2)

    fig.update_layout(barmode="overlay")
    fig.update_traces(opacity=0.7, row=1, col=1)
    fig.update_layout(title_text="Histogram and Boxplot with classes for " + x)
    fig.show()

  elif type=='lineplot':
    fig = px.line(data, x=x, y=y, color=cl)
    fig.show()

  else:
    raise NotImplementedError('Not implemented!')

  return 0

def iqr_outliers(df: pd.DataFrame, col: str):
  """
  Function for managin outliers in given column from dataframe using IQR method. 
  In descriptive statistics, the interquartile range (IQR) is a measure of statistical dispersion, which is the spread of the data. 
  It is defined as the difference between the 75th and 25th percentiles of the data.

  Args:
    df - dataset to remove outliers

    col - column for removing outliers
    
  Returns:
    Whole dataframe with specified column cleaned from outliers using IQR method

  """
  Q1 = df[col].quantile(0.25)
  Q3 = df[col].quantile(0.75)
  IQR = Q3 - Q1
  df.loc[df[col] <= Q1 - 1.5*IQR, col] = df[col].quantile(0.1)
  df.loc[df[col] >= Q3 + 1.5*IQR, col] = df[col].quantile(0.9)

  return df

def balance_data(X, y, algorithm, cat_features, strategy=1, seed = 123, k_neighbors = 5, n_neighbors=3):
  """
  Function for balancing dataset using specified algorithm and sampling strategy. 

  Args:
    X - matrix of features

    y - vector of target variable observations

    algorithm - algorithm for balancing dataset
      values: {'SMOTE', 'SMOTETOMEK', 'SMOTEENN'}

    strategy - sampling information to resample the data set

    Optional **kwargs:
      k_neighbors - number of neighbours for SMOTE algorithm
      n_neighbors - number of neighbours for ADASYN and ENN algorithms

  Returns:
    Whole balanced dataframe 
  """

  ft = []
  for feat in cat_features:
    ft.append(X.columns.get_loc(feat))

  #sm = SMOTE(sampling_strategy = strategy, k_neighbors = k_neighbors)
  sm = SMOTENC(random_state=seed, categorical_features=ft, sampling_strategy = strategy, k_neighbors=k_neighbors)

  if algorithm == 'SMOTE':
    bs = sm

  elif algorithm == 'SMOTETOMEK':
    
    bs = SMOTETomek(random_state = seed, sampling_strategy = strategy, smote= sm)

  elif algorithm == 'SMOTEENN':
    enn = EditedNearestNeighbours(n_neighbors = n_neighbors)

    bs = SMOTEENN(sampling_strategy = strategy, smote=sm, enn = enn)
  
  else:
    raise NotImplementedError('Not implemented!')

  X_bs,  y_bs = bs.fit_resample(X, y)

  df_bs = X_bs.copy()
  df_bs['stroke'] = y_bs

  X_bs = X_bs.drop(columns='index', axis=1)

  print("Liczba wszystkich obserwacji: ", len(y_bs))
  print("\nLiczba obserwacji z klasy mniejszościowej: ", len(y_bs[y_bs==1]))
  print("\nProcent klasy mniejszościowej do całości po zbalansowaniu: ", len(y_bs[y_bs == 1])/len(y_bs))
  
  return df_bs, X_bs, y_bs

#Stwórz klase z doborem parametru thresholda p

class LogisticRegressionWithThreshold(LogisticRegression):
    def predict(self, X, threshold=None):
        if threshold == None: # If no threshold passed in, simply call the base class predict, effectively threshold=0.5
            return LogisticRegression.predict(self, X)
        else:
            y_scores = LogisticRegression.predict_proba(self, X)[:, 1]
            y_pred_with_threshold = (y_scores >= threshold).astype(int)

            return y_pred_with_threshold

class SupportVectorMachineWithThreshold(SVC):
    def predict(self, X, threshold=None):
        if threshold == None: # If no threshold passed in, simply call the base class predict, effectively threshold=0.5
            return SVC.predict(self, X)
        else:
            y_scores = SVC.predict_proba(self, X)[:, 1]
            y_pred_with_threshold = (y_scores >= threshold).astype(int)

            return y_pred_with_threshold

class RandomForestClassifierWithThreshold(RandomForestClassifier):
    def predict(self, X, threshold=None):
        if threshold == None: # If no threshold passed in, simply call the base class predict, effectively threshold=0.5
            return RandomForestClassifier.predict(self, X)
        else:
            y_scores = RandomForestClassifier.predict_proba(self, X)[:, 1]
            y_pred_with_threshold = (y_scores >= threshold).astype(int)

            return y_pred_with_threshold

def preprocess_data(X_train, X_test, columns_to_standardize, cat_features, scaler=True, ohe=True, test_only=False, sc=None):
  """
  Function preprocesses data by standardizing numerical features and one hot encoding categorical features.

  Args:
    X_train - train dataset of features

    X_test - test dataset of features

    columns_to_standardize - columns to standardize

    cat_features - column for OHE

    scaler - boolean value specifies if standardizing should be performed

    ohe - boolean value specifies if OHE should be performed

    test_only - flag if only test set should be preprocessed

    sc - scaler

  Returns
    X 
      Preprocessed train data

    X_t
      Preprocessed test data
  """
  def ohe(df: pd.DataFrame, cat_features):
    df_ohe = df[cat_features].astype('str')
    dummies = pd.get_dummies(df_ohe)

    df = df.drop(columns=cat_features, axis=1)
    df = pd.merge(df, dummies, left_index=True, right_index=True)
    
    return df
  
  if test_only==True:
    X_t = X_test.copy()

    if scaler == True:
      X_t[columns_to_standardize] = sc.transform(X_t[columns_to_standardize])
  
    if ohe == True:
      X_t = ohe(X_t, cat_features)

    return X_t

  else:

    X = X_train.copy()
    X_t = X_test.copy()

    sc = preprocessing.StandardScaler()
    
    if scaler == True:
      X[columns_to_standardize] = sc.fit_transform(X_train[columns_to_standardize])
      X_t[columns_to_standardize] = sc.transform(X_t[columns_to_standardize])
    
    if ohe == True:
      X = ohe(X, cat_features)
      X_t = ohe(X_t, cat_features)

    return X, X_t, sc

def fit_pred_score(X_train, y_train, X_test, y_test, model, model_name, dataset_name, scaler=True, ohe = False, columns_to_standardize = None, cat_features = None, visualize_train = True, visualize_test = True, show_coefficients = False, threshold = None):
  """
  Function for training and testing given classifier. 
  It plots classification reports for train and test data and might display coefficients of predictors. 
  Returns report in dataframe format for given model and dataset with metrics from classification report.

  Args:

    X_train, y_train, X_test, y_test - training and testing data with associated labels

    model - model to train and test

    model_name - name of the model for the report

    dataset_name - dataset name used for the report

    scaler - if features should be scaled

    columns_to_standardize - columns which should be scaled if scaler = True

    visualize_train - if classification report for training data should be visualized

    visualize_test - if classification report for testing data should be visualized

    show_coefficients - if coefficients of features should be displayed

    threshold - value of 'p' for Logistic Regression model

  Returns:

    report - dataframe with metrics from classification report with associated dataset and model name
  """

  X, X_t, standard_scaler = preprocess_data(X_train, X_test, columns_to_standardize = columns_to_standardize, cat_features = cat_features, scaler = scaler, ohe = ohe)

  model.fit(X, y_train)

  if threshold == None:
    y_pred_train = model.predict(X)
    y_pred = model.predict(X_t)
  else:
    y_pred_train = model.predict(X, threshold = threshold)
    y_pred = model.predict(X_t, threshold = threshold)

  classes=['No stroke', 'Stroke']

  #plot classification report for training data
  if visualize_train == True:
    visualizer = ClassificationReport(
          model, classes=classes)
    visualizer.fit(X, y_train)
    visualizer.score(X, y_train)
    visualizer.show()

  if visualize_test == True:
    #plot classification report for test data
    classifier.classification_report(model, X, y_train, X_t, y_test, classes=classes)

  if show_coefficients == True:
    coef_dict = {}
    for coef, feat in zip(model.coef_[0,:],X.columns):
      print(feat, " = ", coef)

  report = pd.DataFrame()
  for t, y, y_approx in zip(['train','test'], [y_train, y_test], [y_pred_train, y_pred]):
    results = pd.DataFrame(metrics.classification_report(y, y_approx, target_names = classes, output_dict=True))
    results['metric'] = results.index
    results.reset_index(inplace=True, drop=True)
    results['dataset_name'] = dataset_name
    results['model_name'] = model_name
    results['split'] = t
    report = pd.concat([report, results])

  
  
  return report

def transform_report_to_plot(dt : pd.DataFrame, x : str, cl = 'metric', train_test = False):
  """
  Function transforms report dataset for plotting in a way that for F1-score macro avg value is taken, for precision the values for no stroke and for recall the values for stroke. 
  
  Args:
    dt - filtered dataset

    x - axis

    cl - class to display

    train_test - flag which specifies if data should be transformed to differences between train and test

  Returns:
    df - transformed report dataframe
  """
  if train_test == True:
    prefix = 'train_test_diff'
  else:
    prefix = ''
  #cols_to_select
  f1_s = dt.loc[(dt.metric == 'f1-score'),[prefix+'macro avg', cl, x]]
  precision_s = dt.loc[(report.metric == 'precision'),[prefix+'Stroke', cl, x]]
  precision_s[cl] = "Precision (Stroke)"
  recall_s = dt.loc[(dt.metric == 'recall'),[prefix+'Stroke', cl, x]]
  recall_s[cl] = "Recall (Stroke)"
  precision_ns = dt.loc[(report.metric == 'precision'),[prefix+'No stroke', cl, x]]
  precision_ns[cl] = "Precision (No stroke)"
  recall_ns = dt.loc[(dt.metric == 'recall'),[prefix+'No stroke', cl, x]]
  recall_ns[cl] = "Recall (No stroke)"

  columns = ['value', cl, '']
  for d in [f1_s, precision_s, recall_s, precision_ns, recall_ns]:
    d.columns = columns

  df = pd.concat([f1_s, precision_s, recall_s, precision_ns, recall_ns])
  return df

def prepare_diff_dataset(df : pd.DataFrame(), by: str, cols_for_diff = ['No stroke',	'Stroke',	'accuracy',	'macro avg',	'weighted avg']):
  """
  Takes the report dataset and returns dataset with calculated differences in metric between test and train by specified grouping.

  Args:

    df - dataset 

    by - specified grouping

    cols_for_diff - cols to calculate differences

  Returns:
    df_diff
      Dataset with calculated differences
  """
  df_diff = df.copy()

  for col in cols_for_diff:
    df_diff['train_test_diff'+col] = None
    for ds in df_diff[by].unique():
      for m in df_diff.metric.unique():
        filter = (df_diff[by] == ds) & (df_diff.metric == m)
        df_diff.loc[filter, 'train_test_diff'+col] = np.abs(df_diff.loc[filter, col].diff())

  return df_diff

Dane

powrót do spisu treści

url = 'https://github.com/maju116/Statistics_II_GUT/raw/main/PROJECT/healthcare-dataset-stroke-data.csv'
df = pd.read_csv(url)

df.head(5)

	id	gender	age	hypertension	heart_disease	ever_married	work_type	Residence_type	avg_glucose_level	bmi	smoking_status	stroke
0	9046	Male	67.0	0	1	Yes	Private	Urban	228.69	36.6	formerly smoked	1
1	51676	Female	61.0	0	0	Yes	Self-employed	Rural	202.21	NaN	never smoked	1
2	31112	Male	80.0	0	1	Yes	Private	Rural	105.92	32.5	never smoked	1
3	60182	Female	49.0	0	0	Yes	Private	Urban	171.23	34.4	smokes	1
4	1665	Female	79.0	1	0	Yes	Self-employed	Rural	174.12	24.0	never smoked	1

df.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 5110 entries, 0 to 5109
Data columns (total 12 columns):
 #   Column             Non-Null Count  Dtype  
---  ------             --------------  -----  
 0   id                 5110 non-null   int64  
 1   gender             5110 non-null   object 
 2   age                5110 non-null   float64
 3   hypertension       5110 non-null   int64  
 4   heart_disease      5110 non-null   int64  
 5   ever_married       5110 non-null   object 
 6   work_type          5110 non-null   object 
 7   Residence_type     5110 non-null   object 
 8   avg_glucose_level  5110 non-null   float64
 9   bmi                4909 non-null   float64
 10  smoking_status     5110 non-null   object 
 11  stroke             5110 non-null   int64  
dtypes: float64(3), int64(4), object(5)
memory usage: 479.2+ KB

df.describe().transpose()

	count	mean	std	min	25%	50%	75%	max
id	5110.0	36517.829354	21161.721625	67.00	17741.250	36932.000	54682.00	72940.00
age	5110.0	43.226614	22.612647	0.08	25.000	45.000	61.00	82.00
hypertension	5110.0	0.097456	0.296607	0.00	0.000	0.000	0.00	1.00
heart_disease	5110.0	0.054012	0.226063	0.00	0.000	0.000	0.00	1.00
avg_glucose_level	5110.0	106.147677	45.283560	55.12	77.245	91.885	114.09	271.74
bmi	4909.0	28.893237	7.854067	10.30	23.500	28.100	33.10	97.60
stroke	5110.0	0.048728	0.215320	0.00	0.000	0.000	0.00	1.00

Mamy łącznie 12 kolumn po 5110 obserwacji z czego jedna z nich stanowi ('stroke') zmienną celu i będziemy starać się ja zamodelować. Jest to zmienna binarna informująca czy dany pacjent miał udar mózgu czy też nie zatem jest to typowe zadanie klasyfikujące. Oprócz tego jest pozostałych 11 zmiennych:

• zmienna id typu integer przypisująca konkretnego pacjenta nie wnosi raczej żadnej wartości informacyjnej ale pozwoli nam sprawdzić czy jakiś pacjent jest dwukrotnie w historii co mogłoby wskazywać na podwyższone ryzyko udaru.
• zmienna gender typu object jest zmienną binarną wskazująca na płeć pacjenta, z pewnością będzie trzeba ją zakodować. Spodziewa się, iż będzie dawała istotną informacje dla zmiennej celu.
• zmienna age typu float jest zmienną numeryczną wskazująca na wiek pacjenta i także można się spodziewać że będzie istotna w modelowaniu. Mamy tutaj zakres wartości od 25 do 82 bez raczej mocnych ogonów.
• zmienna hypertension typu integer jest zmienną binarną informująca o nadciśnieniu pacjenta.
• zmienna heart_disease typu integer jest zmienną binarną informująca o występowaniu chorób serca.
• zmienna ever_married typu object jest zmienną kategoryczną informująca o tym czy pacjent kiedykolwiek był w związku małżeńskim.
• zmienna work_type typu object jest zmienną kategoryczną informująca o typie zatrudnieniu.
• zmienna residence_type typu object jest zmienną kategoryczną informująca o strukturze miejsca zamieszkania.
• zmienna avg_glucose_level typu float jest zmienną numeryczną informująca o przeciętnyn poziomie glukozy we krwi.
• zmienna bmi typu float jest zmienną numeryczną informująca o wskaźniku BMI pacjenta. Posiada 201 braków danych.
• zmienna smoking_status typu object jest zmienną kategoryczną informująca o tym czy i jak często palił pacjent.
• zmienna stroke typu int jest zmienną binarną informująca czy pacjent doznał udaru mózgu i jest to nasza zmienną celu.

Widzimy, że mamy 211 braków danych dla zmiennej bmi, którymi zajmiemy się w dalszej części projektu.

EDA

powrót do spisu treści

stroke

powrót do spisu treści

df.stroke.value_counts()

0    4861
1     249
Name: stroke, dtype: int64

plotting(df, 'stroke', type = 'hist')

0

Widzimy, że nasza zmienna celu jest niezbalansowana gdzie klasą większościową są pacjenci bez udaru mózgu. Musimy mieć to na uwadze i w przypadku pozostawienia takiego stanu rzeczy nie możemy skorzystać z metryki accuracy, a powinniśmy skorzystać z innej - przykładowo średniej harmonicznej między precision a recall czyli f1-score. Niezbalansowany zbiór danych jest jednak dość problemtyczny, gdyż klasyfikator ma tendencję skupiania się na predykcji klasy większościowej.

Istnieją sposoby na zbalansowanie zbioru danych, które są częścią pre-processingu czyli wstępnego przetwarzania danych przed modelowaniem. Do jednych z bardziej popularnych metod należą:

• undersampling - pozostawia wszystkie obserwacje z klasy mniejszościowej i losowo eliminuje obiekty z klasy większościowej
• oversampling - pozostawia wszystkie obserwacje z klasy większościowej i losowo replikuje elementy z klasy mniejszościowej

W przypadku undersamplingu możemy pozbyć się istotnych informacyjnie obserwacji. Dodatkowo nasz zbiór danych jest niewielki, więc nie chcielibyśmy się ograniczyć do zaledwie 249 obserwacji. Natomiast oversampling może doprowadzić do nadmiernego dopasowania modelu do naszych danych. Istnieje też podejście mieszane wykorzystujące zarówno undersampling, jak też oversampling ale ciężko jest ustalić optymalną proporcję między podejściami. Spróbujemy sprawdzić empirycznie najlepszą metodę i się do niej ograniczyć oczywiście opierając się nie o zwykłe próbkowanie, a dostępne i znane algorytmu balansowania zbiorów danych.

id

powrót do spisu treści

print("Liczba zduplikowanych wartości ID pacjentów: ", len(df[df.id.duplicated()]))

Liczba zduplikowanych wartości ID pacjentów:  0

W naszym zbiorze danych każdy pacjent jest wpisany unikalnie, tak więc zmienna ID raczej do niczego więcej nam się nie przyda i będziemy mogli ją pozostawić w dalszych etapach.

gender

powrót do spisu treści

plotting(data = df, x='gender',  type='class_hist')

0

df.gender.value_counts()

Female    2994
Male      2115
Other        1
Name: gender, dtype: int64

W zmiennej gender występują 3 wartości (male, female, other). Wartość 'Other' występuje tylko jeden raz, więc można go potraktować jako brak danych. W zależności od wartości pozostałych zmiennej dla tej obserwacji albo zastąpimy ją wartością male lub female (jako tą która najczęściej występuje) albo usuniemy.

age

powrót do spisu treści

plotting(data = df, x='age',  type='hist_box')

0

plotting(data = df, x='age',  type='class_hist_box')

0

Możemy zaobserwować, że udar mózgu w naszym zbiorze mają głównie osoby w podeszłym wieku (60-80 lat). Natomiast osoby bez udaru rozkładają się dość równomiernie, mediana wynosi 43 lata, gdzie dla osób z udarem 71 lat. Dodatkowo występują zmienne odstające w grupie ludzi z udarem. Sam rozkład zmiennej zbliżony jest do rozkładu jednostajnego z "grubym" oszacowaniem.

hypertension

powrót do spisu treści

plotting(data = df, x='hypertension',  type='class_hist')

0

Można zauważyć, że znaczna większość pacjentów nie ma nadciśnienia. Co więcej, mimo braku zbalansowanego zbioru danych, więcej pacjentów z udarem mózgu nie miało przy tym nadciśnienia.

heart_disease

powrót do spisu treści

plotting(data = df, x='heart_disease',  type='class_hist')

0

Podobna sytuacja jak dla zmiennej hypertension. Znaczna większość pacjentow nie ma chorób serca oraz mimo braku zbalansowania, większość osób z udarem także nie ma chorób serca.

ever_married

powrót do spisu treści

plotting(data = df, x='ever_married',  type='class_hist')

0

Zmienna ever_married jest zmienną kategoryczną, którą będzie trzeba zakodować w wartości (0,1). Większa część pacjentów była kiedykolwiek żonata oraz to w tej grupie udar mózgu występuje w większej ilości.

work_type

powrót do spisu treści

plotting(data = df, x='work_type',  type='class_hist')

0

Zmienna work_type określająca rodzaj wykonywanej pracy jest także zmienną kategoryczną, którą zakodujemy w wartości numeryczne. Najwięcej pacjentów jest z kategorii 'private', następnie równiemiernie ilościowo rozłożeni między 'children', 'self_employed', 'govt_job'. Znikoma ilość nigdy nie pracowała. Największa ilość udarów znajduje się wśród sektora prywatnego oraz samo-zatrunionego. Niewielka część w sektorze publicznym. Pojedyncze wartości występują też w grupie dzieci, ale z poprzednich wniosków najprawdopodobniej są to obserwacje odstające. Można także spróbować zkoszykować tą zmienną i rozdzielić na sektor obciążony i nieobciążony względem ryzyka udaru mózgu.

residence_type

powrót do spisu treści

plotting(data = df, x='Residence_type',  type='class_hist')

0

Zmienna określa obszar zamieszkania z rozdziałem na wiejski i miejski. Widać, że podział jest równomierny oraz w każdym z podziałów jest zbliżona proporcja osób po udarze do osób bez udaru mózgu, więc zmienna ta niesie ze sobą najpewniej niewielką ilość informacji dla potencjalnego modelu.

avg_glucose_level

powrót do spisu treści

plotting(data = df, x='avg_glucose_level',  type='hist_box')

0

plotting(data = df, x='avg_glucose_level',  type='class_hist_box')

0

Możemy zaobserwować, że najwięcej jest pacjentów ze średnim poziomem glukozy w przedziale 77.24 - I kwantyl, 114.09 - II kwantyl, a więc z poziomem prawidłowym lub stanem przedcukrzycowym. Wniosek ten jednak jest uogólniony, gdyż zakres ten jest nieco inny w zależności od pory czy też sposobu wykonywanego badania, a także wieku badanej osoby. Osoby z wyższym średnim poziomem glukozy we krwi są outlierami dla grupy bez udaru mózgu, natomiast dla pozostałej części osób mieszczą się one w IV kwantylu. Więc usunięcie tych obserwacji, lub też zmiana ich wartości spowodowałaby utratę być może znaczącej informacji. Zamiast tego spróbujemy wykorzystać tą zależność tworząc nową zmienną kategoryczną informującą o średnim zakresie glukozy (poziom prawidłowy, stan przedcukrzycowy, cukrzyca) oraz imputując nowe wartości, jak też transformując całą zmienną.

bmi

powrót do spisu treści

plotting(data = df, x='bmi',  type='hist_box')

0

plotting(data = df, x='bmi',  type='class_hist_box')

0

Rozkład zmiennej bmi przypomina rozkład normalny z cięższym prawym ogonem. Występują też outliery dla bmi > 47.5. Większość pacjentów zawiera się w przedziale 23.5 - 33.1 czyli jest z wagą prawidłową lub nadwagą I stopnia. O ile sama zmienna może nie nieść zbyt dużej informacji z uwagi, że rozkład wartości jest zbliżony dla osób z udarem, jak i bez udaru to niemniej możliwe są pewne istotne interakcje między pozostałymi zmiennymi jak chociażby średnim poziomem glukozy we krwi. W kontekście danych odstających również można stworzyć zmienną kategoryczną zgodnie z kategoriami wskaźników BMI.

Dodatkowo zmienna BMI zawiera braki danych, którymi należy się zająć. Istnieje na to wiele metod m.in.

• usunięcie obserwacji
• wypełnienie ich statystyką opisową typu średnia, mediana
• wykorzystanie algorytmów ML takich jak KNN

Z racji, iż braków danych jest niewiele (zaledwie 4%) to wypełnimi ich statystyką opisową jako działaniem prostym, ale i skutecznym. Jak wspomniano, rozkład zmiennej jest zbliżony do normalnego więc nie ma większego znaczenia wybór między medianą, a średnią ale zważywszy na lekki ogon prawostronny wybrana zostanie mediana. Dla obserwacji odstających posłużymy się metodą IQR.

smoking_status

powrót do spisu treści

plotting(data = df, x='smoking_status',  type='class_hist_box')

0

Ciężko wyciągnać tu szersze wnioski zważywszy na dość liczną grupę 'Unknown', Najliczniejsza grupa pacjentów nigdy nie paliła papierosów zakładając, że w grupie 'Unknown' wartości rozkładałyby się w miare równomiernie. Widać rozróżnienei dla osób z udarem, że w większości są to osoby niepalące lub palące sporadycznie. Sprawdźmy wiek osób dla grupy Unknown, być może w naturalny sposób wydzieli nam się grupa osób niepełnoletnich.

plotting(data=df[(df.age <= 16)], x='smoking_status', type='class_hist')

0

Widzimy, że dla osób niepełnoletnich wartość smoking_status jest w większości równa 'Unknown'. Należy coś z tą grupą osób zrobić, jednakże nie można po prostu usunąć tych obserwacji, gdyż udar mózgu dotyczyć może także dzieci. Z jednej strony można zaklasyfikować te osoby do grupy 'Never smoked', jednak jest wysoce prawdopodobne, że w rzeczywistości część osób może 'popalać'. Stworzymy więc nową kategorię 'Underage' dla tych pacjentów.

Inżynieria cech

powrót do spisu treści

W inżynierii cech będziemy tworzyli nowe zmienne z już istniejących, zajmowali się brakiem danych czy też odstającymi obserwacjami. Do tego ostatniego skorzystamy z metody rozstępu międzyćwiartylkowego IQR. Jest to różnica między trzecim a pierwszym kwartylem. Ponieważ pomiędzy tymi kwartylami znajduje się z definicji 50% wszystkich obserwacji (położonych centralnie w rozkładzie), dlatego im większa szerokość rozstępu ćwiartkowego, tym większe zróżnicowanie cechy. Dodatkowo na końcu dokonamy balansowania naszego zbioru z wykorzystaniem 4 technik: SMOTE, SMOTE+Tomek, SMOTE ENN.

#głęboka kopia danych

df_clear = df.copy()

bmi

powrót do spisu treści

• imputacja braków danych przy użyciu mediany
• stworzenie nowej zmiennej kategorycznej na podstawie istniejącej zmiennej numerycznej oraz zakodowanie jej. Do kategoryzacji posłużymy się kategoriami zaproponowanymi przez 'Centers for Disease Control and Prevention' (https://www.cdc.gov/healthyweight/assessing/bmi/adult_bmi/index.html).
• dodanie zmiennej numerycznej bez danych odstających poprzez wykrycie ich metodą IQR oraz ograniczenie ich wartości na 0.9 kwantyl oraz 0.1 kwantyl
• log-transformacja zmiennej oraz wykrycie ich metodą IQR oraz ograniczenie ich wartości na 0.9 kwantyl oraz 0.1 kwantyl

#imputacja mediany w miejsca NA
df_clear['bmi_no_nan'] = df_clear.bmi.fillna(df_clear.bmi.median())

#zmienna numeryczna bez outlierów
df_clear['bmi_no_outliers'] = df_clear.bmi_no_nan
df_clear = iqr_outliers(df_clear, 'bmi_no_outliers')

#log-transformacja bez outlierów
df_clear['log_bmi_no_outliers'] = np.log(df_clear.bmi_no_nan)
df_clear = iqr_outliers(df_clear, 'log_bmi_no_outliers')

#nowa zmienna kategoryczna zgodnie ze wskaźnikami BMI
df_clear.loc[(df_clear.bmi_no_outliers < 18.5), 'bmi_categorical'] = 'underweight'
df_clear.loc[(df_clear.bmi_no_outliers >= 18.5) & (df_clear.bmi_no_outliers < 24.9), 'bmi_categorical'] = 'normal'
df_clear.loc[(df_clear.bmi_no_outliers >= 24.9) & (df_clear.bmi_no_outliers < 29.9), 'bmi_categorical'] = 'overweight'
df_clear.loc[(df_clear.bmi_no_outliers >= 29.9) & (df_clear.bmi_no_outliers < 34.9), 'bmi_categorical'] = 'obese'
df_clear.loc[(df_clear.bmi_no_outliers >= 34.9), 'bmi_categorical'] = 'extremely obese'

#kodowanie medianą przedziałów
df_clear.loc[df_clear.bmi_categorical == 'underweight', 'bmi_categorical_encoded'] = int(df_clear.loc[df_clear.bmi_categorical == 'underweight','bmi_no_outliers'].median())
df_clear.loc[df_clear.bmi_categorical == 'normal', 'bmi_categorical_encoded'] = int(df_clear.loc[df_clear.bmi_categorical == 'normal','bmi_no_outliers'].median())
df_clear.loc[df_clear.bmi_categorical == 'overweight', 'bmi_categorical_encoded'] = int(df_clear.loc[df_clear.bmi_categorical == 'overweight','bmi_no_outliers'].median())
df_clear.loc[df_clear.bmi_categorical == 'obese', 'bmi_categorical_encoded'] = int(df_clear.loc[df_clear.bmi_categorical == 'obese','bmi_no_outliers'].median())
df_clear.loc[df_clear.bmi_categorical == 'extremely obese', 'bmi_categorical_encoded'] = int(df_clear.loc[df_clear.bmi_categorical == 'extremely obese','bmi_no_outliers'].median())
df_clear.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 5110 entries, 0 to 5109
Data columns (total 17 columns):
 #   Column                   Non-Null Count  Dtype  
---  ------                   --------------  -----  
 0   id                       5110 non-null   int64  
 1   gender                   5110 non-null   object 
 2   age                      5110 non-null   float64
 3   hypertension             5110 non-null   int64  
 4   heart_disease            5110 non-null   int64  
 5   ever_married             5110 non-null   object 
 6   work_type                5110 non-null   object 
 7   Residence_type           5110 non-null   object 
 8   avg_glucose_level        5110 non-null   float64
 9   bmi                      4909 non-null   float64
 10  smoking_status           5110 non-null   object 
 11  stroke                   5110 non-null   int64  
 12  bmi_no_nan               5110 non-null   float64
 13  bmi_no_outliers          5110 non-null   float64
 14  log_bmi_no_outliers      5110 non-null   float64
 15  bmi_categorical          5110 non-null   object 
 16  bmi_categorical_encoded  5110 non-null   float64
dtypes: float64(7), int64(4), object(6)
memory usage: 678.8+ KB

for col, types in zip(['bmi_no_nan', 'bmi_categorical', 'bmi_no_outliers', 'log_bmi_no_outliers'], ['class_hist_box','class_hist','class_hist_box','class_hist_box']):
  plotting(data = df_clear, x=col, type=types)

Możemy zaobserwować, że wypełnienie braków danych spowodowało peak w punkcie mediany co jest naturalnym efektem zastosowanej metody.

Dodatkowo, najbardziej liczna jest grupa osób z nadwagą i to tam również najliczniej występują osoby z udarem mózgu. W przypadku usunięcia outlierów mozemy zaobserwować pozbycie się ogona prawostronnego, natomiast zastosowanie log-transformacji wydaje się nie mieć znacznego wpływu poza przeskalowaniem co będzie zbędne z racji iż w późniejszym etapie będziemy dokonywać standaryzacji zmiennych, dlatego odrzucimy tą zmienną.

Co więcej, możemy zauważyć, że przypisanie granicznych wartości kwantylowych nie zlikwidowało danych odstających w grupie osób z udarem. W tej sytuacji możemy usunąć obserwacje specjalnie dla tej grupy, ale nie są to pojedyncze przypadki a i tak jest ona mało liczna. W takim razie, z racji dalszego lekkiego zaszumienia zmiennej przez te obserwacje, pozostaniemy przy zmiennej kategorycznej koszykującej do odpowiedniej grupy. Jest to dość powszechne rozwiązanie dla BMI, gdyż to powszechny indeks wagi ciała.

gender

powrót do spisu treści

• z racji, iż jest to tylko jedna obserwacja i ciężko jest ją odpowiednio podmienić (pozostałe zmienne są mocno rozrzucone, więc nie niosą wystarczającej informacji do imputacji) zostanie usunięta
• kodowanie zmiennej

df_clear = df_clear.loc[df_clear.gender != 'Other',:]
df_clear.loc[df_clear.gender == 'Female', 'is_female'] = 1
df_clear.loc[df_clear.gender == 'Male', 'is_female'] = 0
df_clear.is_female = df_clear.is_female.astype('int')

df_clear.info()

<class 'pandas.core.frame.DataFrame'>
Int64Index: 5109 entries, 0 to 5109
Data columns (total 18 columns):
 #   Column                   Non-Null Count  Dtype  
---  ------                   --------------  -----  
 0   id                       5109 non-null   int64  
 1   gender                   5109 non-null   object 
 2   age                      5109 non-null   float64
 3   hypertension             5109 non-null   int64  
 4   heart_disease            5109 non-null   int64  
 5   ever_married             5109 non-null   object 
 6   work_type                5109 non-null   object 
 7   Residence_type           5109 non-null   object 
 8   avg_glucose_level        5109 non-null   float64
 9   bmi                      4908 non-null   float64
 10  smoking_status           5109 non-null   object 
 11  stroke                   5109 non-null   int64  
 12  bmi_no_nan               5109 non-null   float64
 13  bmi_no_outliers          5109 non-null   float64
 14  log_bmi_no_outliers      5109 non-null   float64
 15  bmi_categorical          5109 non-null   object 
 16  bmi_categorical_encoded  5109 non-null   float64
 17  is_female                5109 non-null   int32  
dtypes: float64(7), int32(1), int64(4), object(6)
memory usage: 738.4+ KB

age

powrót do spisu treści

• pozbycie się outlierów dla klasy 'stroke', gdyż widać zróżnicowanie dla tej grupy i pozostawienie ich mogłoby wprowadzić zbędny szum

#tablica z pacjentami z udarem
df_stroke = df_clear[df_clear.stroke == 1].copy()
#pozbycie sie outlierow
df_stroke = iqr_outliers(df=df_stroke, col='age')
#ponowna konkatenacja zbioru
df_clear = pd.concat([df_clear[df_clear.stroke==0], df_stroke])
#plot
plotting(data=df_clear, x='age', type='class_hist_box')
#podsumowanie
df_clear.info()

<class 'pandas.core.frame.DataFrame'>
Int64Index: 5109 entries, 249 to 248
Data columns (total 18 columns):
 #   Column                   Non-Null Count  Dtype  
---  ------                   --------------  -----  
 0   id                       5109 non-null   int64  
 1   gender                   5109 non-null   object 
 2   age                      5109 non-null   float64
 3   hypertension             5109 non-null   int64  
 4   heart_disease            5109 non-null   int64  
 5   ever_married             5109 non-null   object 
 6   work_type                5109 non-null   object 
 7   Residence_type           5109 non-null   object 
 8   avg_glucose_level        5109 non-null   float64
 9   bmi                      4908 non-null   float64
 10  smoking_status           5109 non-null   object 
 11  stroke                   5109 non-null   int64  
 12  bmi_no_nan               5109 non-null   float64
 13  bmi_no_outliers          5109 non-null   float64
 14  log_bmi_no_outliers      5109 non-null   float64
 15  bmi_categorical          5109 non-null   object 
 16  bmi_categorical_encoded  5109 non-null   float64
 17  is_female                5109 non-null   int32  
dtypes: float64(7), int32(1), int64(4), object(6)
memory usage: 738.4+ KB

ever_married

powrót do spisu treści

• zakodowanie zmiennej

df_clear.loc[df_clear.ever_married == 'Yes', 'is_married'] = 1
df_clear.loc[df_clear.ever_married == 'No', 'is_married'] = 0
df_clear.is_married = df_clear.is_married.astype('int')

df_clear.info()

<class 'pandas.core.frame.DataFrame'>
Int64Index: 5109 entries, 249 to 248
Data columns (total 19 columns):
 #   Column                   Non-Null Count  Dtype  
---  ------                   --------------  -----  
 0   id                       5109 non-null   int64  
 1   gender                   5109 non-null   object 
 2   age                      5109 non-null   float64
 3   hypertension             5109 non-null   int64  
 4   heart_disease            5109 non-null   int64  
 5   ever_married             5109 non-null   object 
 6   work_type                5109 non-null   object 
 7   Residence_type           5109 non-null   object 
 8   avg_glucose_level        5109 non-null   float64
 9   bmi                      4908 non-null   float64
 10  smoking_status           5109 non-null   object 
 11  stroke                   5109 non-null   int64  
 12  bmi_no_nan               5109 non-null   float64
 13  bmi_no_outliers          5109 non-null   float64
 14  log_bmi_no_outliers      5109 non-null   float64
 15  bmi_categorical          5109 non-null   object 
 16  bmi_categorical_encoded  5109 non-null   float64
 17  is_female                5109 non-null   int32  
 18  is_married               5109 non-null   int32  
dtypes: float64(7), int32(2), int64(4), object(6)
memory usage: 758.4+ KB

work_type

powrót do spisu treści

• zakodowanie zmiennej

le = preprocessing.LabelEncoder()
le.fit(df_clear.work_type)
df_clear['work_type_encoded'] = le.transform(df_clear.work_type)
df_clear.info()

<class 'pandas.core.frame.DataFrame'>
Int64Index: 5109 entries, 249 to 248
Data columns (total 20 columns):
 #   Column                   Non-Null Count  Dtype  
---  ------                   --------------  -----  
 0   id                       5109 non-null   int64  
 1   gender                   5109 non-null   object 
 2   age                      5109 non-null   float64
 3   hypertension             5109 non-null   int64  
 4   heart_disease            5109 non-null   int64  
 5   ever_married             5109 non-null   object 
 6   work_type                5109 non-null   object 
 7   Residence_type           5109 non-null   object 
 8   avg_glucose_level        5109 non-null   float64
 9   bmi                      4908 non-null   float64
 10  smoking_status           5109 non-null   object 
 11  stroke                   5109 non-null   int64  
 12  bmi_no_nan               5109 non-null   float64
 13  bmi_no_outliers          5109 non-null   float64
 14  log_bmi_no_outliers      5109 non-null   float64
 15  bmi_categorical          5109 non-null   object 
 16  bmi_categorical_encoded  5109 non-null   float64
 17  is_female                5109 non-null   int32  
 18  is_married               5109 non-null   int32  
 19  work_type_encoded        5109 non-null   int32  
dtypes: float64(7), int32(3), int64(4), object(6)
memory usage: 778.3+ KB

residence_type

powrót do spisu treści

• zakodowanie zmiennej

df_clear.loc[df_clear.Residence_type == 'Rural', 'is_rural'] = 1
df_clear.loc[df_clear.Residence_type == 'Urban', 'is_rural'] = 0 
df_clear.is_rural = df_clear.is_rural.astype('int')

df_clear.info()

<class 'pandas.core.frame.DataFrame'>
Int64Index: 5109 entries, 249 to 248
Data columns (total 21 columns):
 #   Column                   Non-Null Count  Dtype  
---  ------                   --------------  -----  
 0   id                       5109 non-null   int64  
 1   gender                   5109 non-null   object 
 2   age                      5109 non-null   float64
 3   hypertension             5109 non-null   int64  
 4   heart_disease            5109 non-null   int64  
 5   ever_married             5109 non-null   object 
 6   work_type                5109 non-null   object 
 7   Residence_type           5109 non-null   object 
 8   avg_glucose_level        5109 non-null   float64
 9   bmi                      4908 non-null   float64
 10  smoking_status           5109 non-null   object 
 11  stroke                   5109 non-null   int64  
 12  bmi_no_nan               5109 non-null   float64
 13  bmi_no_outliers          5109 non-null   float64
 14  log_bmi_no_outliers      5109 non-null   float64
 15  bmi_categorical          5109 non-null   object 
 16  bmi_categorical_encoded  5109 non-null   float64
 17  is_female                5109 non-null   int32  
 18  is_married               5109 non-null   int32  
 19  work_type_encoded        5109 non-null   int32  
 20  is_rural                 5109 non-null   int32  
dtypes: float64(7), int32(4), int64(4), object(6)
memory usage: 798.3+ KB

avg_glucose_level

powrót do spisu treści

• stworzenie nowej zmiennej kategorycznej na podstawie istniejącej zmiennej numerycznej oraz zakodowanie jej
• dodanie zmiennej numerycznej bez danych odstających poprzez wykrycie ich metodą IQR oraz ograniczenie ich wartości na 0.9 kwantyl oraz 0.1 kwantyl
• log-transformacja zmiennej oraz wykrycie ich metodą IQR oraz ograniczenie ich wartości na 0.9 kwantyl oraz 0.1 kwantyl

#nowa zmienna kategoryczna zgodnie ze wskaźnikami BMI
df_clear.loc[(df_clear.avg_glucose_level < 140), 'glucose_categorical'] = 'normal'
df_clear.loc[(df_clear.avg_glucose_level >= 140) & (df_clear.avg_glucose_level <= 199), 'glucose_categorical'] = 'prediabetes'
df_clear.loc[(df_clear.avg_glucose_level > 199), 'glucose_categorical'] = 'diabetes'

#kodowanie
le_gl = preprocessing.LabelEncoder()
le_gl.fit(df_clear.glucose_categorical)
df_clear['glucose_categorical_encoded'] = le_gl.transform(df_clear.glucose_categorical)
df_clear.info()

#zmienna numeryczna bez outlierów
df_clear['glucose_no_outliers'] = df_clear.avg_glucose_level
df_clear = iqr_outliers(df_clear, 'glucose_no_outliers')

#log-transformacja bez outlierów
df_clear['log_glucose_no_outliers'] = np.log(df_clear.avg_glucose_level)
df_clear = iqr_outliers(df_clear, 'log_glucose_no_outliers')

<class 'pandas.core.frame.DataFrame'>
Int64Index: 5109 entries, 249 to 248
Data columns (total 23 columns):
 #   Column                       Non-Null Count  Dtype  
---  ------                       --------------  -----  
 0   id                           5109 non-null   int64  
 1   gender                       5109 non-null   object 
 2   age                          5109 non-null   float64
 3   hypertension                 5109 non-null   int64  
 4   heart_disease                5109 non-null   int64  
 5   ever_married                 5109 non-null   object 
 6   work_type                    5109 non-null   object 
 7   Residence_type               5109 non-null   object 
 8   avg_glucose_level            5109 non-null   float64
 9   bmi                          4908 non-null   float64
 10  smoking_status               5109 non-null   object 
 11  stroke                       5109 non-null   int64  
 12  bmi_no_nan                   5109 non-null   float64
 13  bmi_no_outliers              5109 non-null   float64
 14  log_bmi_no_outliers          5109 non-null   float64
 15  bmi_categorical              5109 non-null   object 
 16  bmi_categorical_encoded      5109 non-null   float64
 17  is_female                    5109 non-null   int32  
 18  is_married                   5109 non-null   int32  
 19  work_type_encoded            5109 non-null   int32  
 20  is_rural                     5109 non-null   int32  
 21  glucose_categorical          5109 non-null   object 
 22  glucose_categorical_encoded  5109 non-null   int32  
dtypes: float64(7), int32(5), int64(4), object(7)
memory usage: 858.2+ KB

for col, types in zip(['glucose_categorical', 'glucose_no_outliers', 'log_glucose_no_outliers'], ['class_hist','class_hist_box','class_hist_box']):
  plotting(data = df_clear, x=col, type=types)

Najbardziej liczna jest grupa osób z prawidłowym kategorią średniego poziomu glukozy we krwi i to tam również najliczniej występują osoby z udarem mózgu (zaliczono tam także I stopień otyłości).

W kontekście kategoryzacji zmiennej otrzymano raczej spadek informacji, gdyż została jedna najliczniejsza grupa osób z udarem oraz pozostałe z podobną licznością. Co więcej, jak wspomniano, ciężko jednoznacznie wpisać w kategorie gdyż nie znamy właściwości pomiaru cukru, jak też dodatkowych chorób pacjentów takich jak cukrzyca.

Widzimy także, że metoda IQR z imputacją zamiast usunięcia danych dalej zachowała pewne dane odstające. Jednakże przy skorzystaniu z log-transformacji, są to nieliczne wartości mocno zbliżone do wartości górnej. Mając to na uwadze, jak też fakt spadku informacji po kategoryzacji ograniczymy się do tej zmiennej.

smoking_status

powrót do spisu treści

• stworzenie nowej grupy 'Underage' dla osób niepełnoletnich
• zakodowanie zmiennej

df_clear['smoking_status_new'] = df_clear.smoking_status
df_clear.loc[df_clear.age <= 21, 'smoking_status_new'] = 'underage'
print("Grupa 'Unknown' stanowi %.2f procent obserwacji całościowego zbioru" % (len(df_clear[df_clear.smoking_status_new == 'underage'])/len(df_clear)))
plotting(data=df_clear, x='smoking_status_new', type='class_hist')
df_clear.info()

Grupa 'Unknown' stanowi 0.21 procent obserwacji całościowego zbioru

<class 'pandas.core.frame.DataFrame'>
Int64Index: 5109 entries, 249 to 248
Data columns (total 26 columns):
 #   Column                       Non-Null Count  Dtype  
---  ------                       --------------  -----  
 0   id                           5109 non-null   int64  
 1   gender                       5109 non-null   object 
 2   age                          5109 non-null   float64
 3   hypertension                 5109 non-null   int64  
 4   heart_disease                5109 non-null   int64  
 5   ever_married                 5109 non-null   object 
 6   work_type                    5109 non-null   object 
 7   Residence_type               5109 non-null   object 
 8   avg_glucose_level            5109 non-null   float64
 9   bmi                          4908 non-null   float64
 10  smoking_status               5109 non-null   object 
 11  stroke                       5109 non-null   int64  
 12  bmi_no_nan                   5109 non-null   float64
 13  bmi_no_outliers              5109 non-null   float64
 14  log_bmi_no_outliers          5109 non-null   float64
 15  bmi_categorical              5109 non-null   object 
 16  bmi_categorical_encoded      5109 non-null   float64
 17  is_female                    5109 non-null   int32  
 18  is_married                   5109 non-null   int32  
 19  work_type_encoded            5109 non-null   int32  
 20  is_rural                     5109 non-null   int32  
 21  glucose_categorical          5109 non-null   object 
 22  glucose_categorical_encoded  5109 non-null   int32  
 23  glucose_no_outliers          5109 non-null   float64
 24  log_glucose_no_outliers      5109 non-null   float64
 25  smoking_status_new           5109 non-null   object 
dtypes: float64(9), int32(5), int64(4), object(8)
memory usage: 977.9+ KB

W dalszym ciągu grupa Unknown jest dość liczna, ale wydobyliśmy za to bardziej zróżnicowany podział i pozostawimy to w ten sposób. Dzięki temu możemy dostrzec, że w prawie każdej grupie mamy równy proporcję osób z udarem, jak też bez udaru. Nie jest to zbyt informacyjna zmienna wyszczególniająca informację na temat pacjentów z udarem. Moglibyśmy tutaj wyodrebnić grupę 'underage', a pozostałe zagregować otrzymując zmienną binarną, jednak będzie ona opierała się tak naprawdę tylko na wieku, a nie na rzeczywistym statusie palenia. Odrzucimy więc tą zmienną w dalszej analizie.

Analiza korelacji

powrót do spisu treści

# exclude categorical data
plotting(data=df_clear.loc[:,~df_clear.columns.isin(
    ['id', 'bmi_categorical_encoded','work_type_encoded','glucose_categorical_encoded'])], x='none', type='heatmap')

0

Widoczna jest słaba korelacja między zmiennymi oraz między zmiennymi, a zmienną celu. Największe wartości widoczne są przy tych samych zmiennych, ale po przeprowadzeniu transformacji lub kodowania co jest rezultatem oczywistym - jedynie zmienna is_married koreluje ze zmienną age ale nie będziemy jednak wykluczać żadnej z nich. Nie występuje więc tutaj problem współliniowości. Ciężko także na tej podstawie dokonać selekcji zmiennych, więc raczej skorzystamy w późniejszym etapie z regularyzacji LASSO.

Polynomial Features

powrót do spisu treści

Bardzo często wejściowe zmienne są interaktywne ze sobą w sposób nieliniowy. Takie interkacje mogą zostać później wykryte i wykorzystane przez model przy etapie selekcji zmiennych. Co więcej, podniesienie niektórych zmiennych do pewnej potęgi może pomóc w wyodrębnieniu istotnych relacji między tą zmienną, a zmienną celu. Przeważnie wystarczy podniesienie do stopnia drugiego (ewentualnie trzeciego), gdyż wielomiany stopnia czwartego i wyższego są nadmiernie elastyczne, mają nieregularne kształty i mogą prowadzić do overfittingu. My wykorzystamy cechy wielomianowe stopnia drugiego.

#get PolynomialFeatures transfomer
poly = preprocessing.PolynomialFeatures(2)

#select only apriori chosen variables
X_poly = df_clear.select_dtypes(['int32','int64','float32','float64'])
X_poly = X_poly.loc[:,~X_poly.columns.isin(['id', 'bmi_no_nan','bmi','avg_glucose_level','glucose_no_outliers', 'log_bmi_no_outliers', 'bmi_no_outliers','glucose_categorical_encoded'])]
print("Ramka wejściowa")
X_poly.info()

print("\nRamka wyjściowa po PolyFeatures transformacji")
df_poly = pd.DataFrame(poly.fit_transform(X_poly), columns = poly.get_feature_names_out(X_poly.columns))
df_poly.info()

Ramka wejściowa
<class 'pandas.core.frame.DataFrame'>
Int64Index: 5109 entries, 249 to 248
Data columns (total 10 columns):
 #   Column                   Non-Null Count  Dtype  
---  ------                   --------------  -----  
 0   age                      5109 non-null   float64
 1   hypertension             5109 non-null   int64  
 2   heart_disease            5109 non-null   int64  
 3   stroke                   5109 non-null   int64  
 4   bmi_categorical_encoded  5109 non-null   float64
 5   is_female                5109 non-null   int32  
 6   is_married               5109 non-null   int32  
 7   work_type_encoded        5109 non-null   int32  
 8   is_rural                 5109 non-null   int32  
 9   log_glucose_no_outliers  5109 non-null   float64
dtypes: float64(3), int32(4), int64(3)
memory usage: 359.2 KB

Ramka wyjściowa po PolyFeatures transformacji
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 5109 entries, 0 to 5108
Data columns (total 66 columns):
 #   Column                                           Non-Null Count  Dtype  
---  ------                                           --------------  -----  
 0   1                                                5109 non-null   float64
 1   age                                              5109 non-null   float64
 2   hypertension                                     5109 non-null   float64
 3   heart_disease                                    5109 non-null   float64
 4   stroke                                           5109 non-null   float64
 5   bmi_categorical_encoded                          5109 non-null   float64
 6   is_female                                        5109 non-null   float64
 7   is_married                                       5109 non-null   float64
 8   work_type_encoded                                5109 non-null   float64
 9   is_rural                                         5109 non-null   float64
 10  log_glucose_no_outliers                          5109 non-null   float64
 11  age^2                                            5109 non-null   float64
 12  age hypertension                                 5109 non-null   float64
 13  age heart_disease                                5109 non-null   float64
 14  age stroke                                       5109 non-null   float64
 15  age bmi_categorical_encoded                      5109 non-null   float64
 16  age is_female                                    5109 non-null   float64
 17  age is_married                                   5109 non-null   float64
 18  age work_type_encoded                            5109 non-null   float64
 19  age is_rural                                     5109 non-null   float64
 20  age log_glucose_no_outliers                      5109 non-null   float64
 21  hypertension^2                                   5109 non-null   float64
 22  hypertension heart_disease                       5109 non-null   float64
 23  hypertension stroke                              5109 non-null   float64
 24  hypertension bmi_categorical_encoded             5109 non-null   float64
 25  hypertension is_female                           5109 non-null   float64
 26  hypertension is_married                          5109 non-null   float64
 27  hypertension work_type_encoded                   5109 non-null   float64
 28  hypertension is_rural                            5109 non-null   float64
 29  hypertension log_glucose_no_outliers             5109 non-null   float64
 30  heart_disease^2                                  5109 non-null   float64
 31  heart_disease stroke                             5109 non-null   float64
 32  heart_disease bmi_categorical_encoded            5109 non-null   float64
 33  heart_disease is_female                          5109 non-null   float64
 34  heart_disease is_married                         5109 non-null   float64
 35  heart_disease work_type_encoded                  5109 non-null   float64
 36  heart_disease is_rural                           5109 non-null   float64
 37  heart_disease log_glucose_no_outliers            5109 non-null   float64
 38  stroke^2                                         5109 non-null   float64
 39  stroke bmi_categorical_encoded                   5109 non-null   float64
 40  stroke is_female                                 5109 non-null   float64
 41  stroke is_married                                5109 non-null   float64
 42  stroke work_type_encoded                         5109 non-null   float64
 43  stroke is_rural                                  5109 non-null   float64
 44  stroke log_glucose_no_outliers                   5109 non-null   float64
 45  bmi_categorical_encoded^2                        5109 non-null   float64
 46  bmi_categorical_encoded is_female                5109 non-null   float64
 47  bmi_categorical_encoded is_married               5109 non-null   float64
 48  bmi_categorical_encoded work_type_encoded        5109 non-null   float64
 49  bmi_categorical_encoded is_rural                 5109 non-null   float64
 50  bmi_categorical_encoded log_glucose_no_outliers  5109 non-null   float64
 51  is_female^2                                      5109 non-null   float64
 52  is_female is_married                             5109 non-null   float64
 53  is_female work_type_encoded                      5109 non-null   float64
 54  is_female is_rural                               5109 non-null   float64
 55  is_female log_glucose_no_outliers                5109 non-null   float64
 56  is_married^2                                     5109 non-null   float64
 57  is_married work_type_encoded                     5109 non-null   float64
 58  is_married is_rural                              5109 non-null   float64
 59  is_married log_glucose_no_outliers               5109 non-null   float64
 60  work_type_encoded^2                              5109 non-null   float64
 61  work_type_encoded is_rural                       5109 non-null   float64
 62  work_type_encoded log_glucose_no_outliers        5109 non-null   float64
 63  is_rural^2                                       5109 non-null   float64
 64  is_rural log_glucose_no_outliers                 5109 non-null   float64
 65  log_glucose_no_outliers^2                        5109 non-null   float64
dtypes: float64(66)
memory usage: 2.6 MB

Jak widzimy uzyskaliśmy wszystkie możliwe kombinacje zmiennych wielomianowych 2-giego stopnia. Wykluczymy jednak na start zmienne, których z pewnością nie wykorzystamy takie jak bias: '1' oraz zmienne binarne podniesione do kwadratu, jak też zmienne z interakcją 'stroke' w celu uniknięcia tzw. data leakage.

col_drop = list(df_poly.filter(regex='stroke'))[1:] + ['1','hypertension^2','heart_disease^2','is_female^2','is_married^2','is_rural^2']
df_poly = df_poly[df_poly.columns.drop(col_drop)]

cat_features = ['is_rural','is_married', 'is_female', 'hypertension', 'heart_disease', 'work_type_encoded']

col_filter = set(list(df_poly.filter(regex='age')) + list(df_poly.filter(regex='log_glucose_no_outliers')) + cat_features + ['bmi_categorical_encoded','stroke'])
df_poly = df_poly[col_filter]
df_poly.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 5109 entries, 0 to 5108
Data columns (total 27 columns):
 #   Column                                           Non-Null Count  Dtype  
---  ------                                           --------------  -----  
 0   heart_disease log_glucose_no_outliers            5109 non-null   float64
 1   is_rural                                         5109 non-null   float64
 2   work_type_encoded                                5109 non-null   float64
 3   age hypertension                                 5109 non-null   float64
 4   is_married log_glucose_no_outliers               5109 non-null   float64
 5   is_female                                        5109 non-null   float64
 6   age bmi_categorical_encoded                      5109 non-null   float64
 7   age is_female                                    5109 non-null   float64
 8   stroke                                           5109 non-null   float64
 9   age^2                                            5109 non-null   float64
 10  age log_glucose_no_outliers                      5109 non-null   float64
 11  age                                              5109 non-null   float64
 12  is_female log_glucose_no_outliers                5109 non-null   float64
 13  is_rural log_glucose_no_outliers                 5109 non-null   float64
 14  bmi_categorical_encoded                          5109 non-null   float64
 15  bmi_categorical_encoded log_glucose_no_outliers  5109 non-null   float64
 16  age is_rural                                     5109 non-null   float64
 17  log_glucose_no_outliers                          5109 non-null   float64
 18  work_type_encoded log_glucose_no_outliers        5109 non-null   float64
 19  log_glucose_no_outliers^2                        5109 non-null   float64
 20  is_married                                       5109 non-null   float64
 21  hypertension                                     5109 non-null   float64
 22  heart_disease                                    5109 non-null   float64
 23  age heart_disease                                5109 non-null   float64
 24  age is_married                                   5109 non-null   float64
 25  age work_type_encoded                            5109 non-null   float64
 26  hypertension log_glucose_no_outliers             5109 non-null   float64
dtypes: float64(27)
memory usage: 1.1 MB

Balansowanie zbioru

powrót do spisu treści

• zbalansowanie zbioru danych przy użyciu algorytmów SMOTE, SMOTE+TOMEK, SMOTE ENN. Przy wykonywaniu balansowania będziemy dążyć do uzyskania proporcji 50/50. W późniejszych etapach ewentualnie powrócimy do najlepszej techniki balansowania i spróbujemy innych proporcji klas.

Ważne jest, aby balansować dane tylko na zbiorze treningowym, gdyż w ogólności zbiór testowy jest zbiorem, który ma odzwierciedlać realne dane - a te zbalansowane nie są.

#wyciagniecie macierzy zmiennych X oraz wektora zmiennej celu Y
X = df_poly.loc[:,~df_poly.columns.isin(['stroke'])]
Y = df_poly.stroke

#train-test split, aby balansować tylko na zbiorze treningowym
X_train, X_test, y_train, y_test = train_test_split(X, Y, test_size=0.33, random_state=123)

#for visualization
X_train['index'] = X_train.reset_index().index

df_imbalanced = X_train.copy()
df_imbalanced['stroke'] = y_train

#plot ukazujacy problem niezbalansowanych danych - klasa pacjentow z udarem jest prawie niewidoczna
plotting(data = df_imbalanced,x='index', y='age', type='scatter')

0

SMOTE

powrót do spisu treści

SMOTE wykorzystuje podejście algorytmiczne KNN tj. dla każdej obserwacji $x$ wybiera $k$ najbliższych sąsiadów i wybiera losowo jednego $x'$ z nich. Następnie różnica między koordynatami $x$ oraz $x'$ jest obliczana i przemnożona losowo przez liczby z przedziału $(0,1)$. Tak stworzona z różnic obserwacja $x''$ jest dodana do zbioru. Geometrycznie odwzorowuje to stworzenie nowej obserwacji na podstawie przesunięcia obserwacji bieżącej do jednego z sąsiadów.

# SMOTE
df_smote, X_smote, y_smote = balance_data(X=X_train, y=y_train, cat_features = cat_features, algorithm = 'SMOTE', strategy=1, k_neighbors = 7)
plotting(data = df_smote, x='index', y='age', type='scatter')

Liczba wszystkich obserwacji:  6530

Liczba obserwacji z klasy mniejszościowej:  3265

Procent klasy mniejszościowej do całości po zbalansowaniu:  0.5

0

SMOTE+Tomek

powrót do spisu treści

SMOTE + Tomek jest techniką hybrydową, która łączy technikę oversamplingu algorytmu SMOTE z techniką undersamplingu algorytmu Tomek. Po zastosowaniu SMOTE'a istnieje możliwość nakładania się klastrów klas co może skutkować w overfittingu. Wówczas użyty zostaje dodatkowo algorytm Tomek, który polega na sparowaniu ze sobą próbek z różnych klas, a następnie obserwacji z klasy większościowej zostają usunięte w celu (niekiedy z obydwu) wyraźniejszej separacji klas.

# SMOTE Tomek
df_smtom, X_smtom, y_smtom = balance_data(X=X_train, y=y_train, cat_features = cat_features, algorithm='SMOTETOMEK', strategy = 1, k_neighbors=7)
plotting(data = df_smtom, x='index', y='age', type='scatter')

Liczba wszystkich obserwacji:  6432

Liczba obserwacji z klasy mniejszościowej:  3216

Procent klasy mniejszościowej do całości po zbalansowaniu:  0.5

0

SMOTE + ENN

powrót do spisu treści

SMOTE + ENN jest techniką hybrydową, która polega na dodatkowym usunięciu pewnej liczby obserwacji z przestrzeni cech. Gdy zbiór zostanie poddany oversamplingu przez SMOTE to następnie zostaje zastosowana technika ENN, która jest techniką undersamplingu. ENN (Edited Nearest Neighbor) polega na odnalezieniu najbliższych K-sąsiadów każdej z obserwacji. Następnie sprawdza czy klasa większościowa tych sąsiadów jest zgodna z klasą danej obserwacji, jeśli nie to zarówno obserwacja jak i sąsiedzi zostają usunięci. Dzięki takiemu podejściu oddzielenie klas od siebie jest wyraźniejsze, aniżeli w pozostałych przypadkach.

# SMOTE ENN
df_smenn, X_smenn, y_smenn = balance_data(X=X_train, y=y_train, cat_features = cat_features, algorithm='SMOTEENN', strategy = 1)
plotting(data = df_smenn, x='index', y='age',type='scatter')

Liczba wszystkich obserwacji:  6320

Liczba obserwacji z klasy mniejszościowej:  3055

Procent klasy mniejszościowej do całości po zbalansowaniu:  0.4833860759493671

0

Widzimy, że każdy z algorytmów balansuje w specyficzny dla siebie sposób generuje sztuczne obserwacje dla klasy mniejszościowej. W algorytmie SMOTE ENN widać dodatkowo zmalenie liczby obserwacji za sprawą wyraźnej techniki undersamplingu. Przeprowadzimy trenowanie modeli na każdym ze zbiorze dla każdego z modeli, aby mieć szerszy widok na rezultaty. Co więcej dla każdego z nich przeprowadzimy także regularyzację LASSO.

Metryki

powrót do spisu treści

W ogólności dla oceny rozwiązań problemów klasyfikacyjnych korzysta się z macierzy pomyłek.

Zestawia ona rzeczywiste wartości klas z wartościami klas zamodelowanymi. W przypadku binarnym można to przedstawić w postaci TP czyli przypadki zaklasyfikowane prawidłowo do klasy prawdziwej, TN czyli przypadki zaklasyfikowane prawidłowo do klasy fałszywej, FP czyli przypadki zaklasyfikowane nieprawidłowo do klasy prawdziwej oraz FN czyli przypadki zaklasyfikowane nieprawidłowo do klasy fałszywej. Wartości FP są tzw. błędem I typu, czyli zaklasyfikowaniem czegoś jako prawdziwe, gdy w rzeczywistości takie nie jest. Z kolei wartości FN są błędem II typu, czyli zaklasyfikowaniem czegoś jako nieprawdziwe, gdy w rzeczywistości prawdziwe jest.

Z takiego zestawienia najpopularniejsze są metryki:

• $Accuracy = \frac{(TP+TN)}{(TP+FP+TN+FN)}$ czyli proporcja poprawnych wyników spośród wszystkich możliwych
• $Precision = \frac{TP}{(TP+FP)}$ czyli precyzja modelu czyli ile z zaklasyfikowanych do klasy prawdziwej rzeczywiście do niej należy
• $Recall = \frac{TP}{(TP+FN)}$ czyli proporcja poprawnych wyników z klasy prawdziwej do wszystkich możliwych do uzyskania
• $F_1 = 2\cdot \frac{precision \cdot recall}{precision + recall}$ czyli średnia harmoniczna pomiędzy metrykami Precision i Recall

W naszym przypadku chcemy się skupić na minimalizacji błędu II typu kosztem błędu typu I. Wynika to z faktu, iż w przypadku popełnienia błędu I typu czyli nieprawdziwej predykcji udaru mózgu mogą zostać wykonane niepotrzebne badania oraz poniesione środki finansowe, natomiast w przypadek błędu II typu niesie poważne konsekwencje zdrowotne a w najgorszym przypadku śmierć.

Nie wystarczy jednak minimalizacja FN, ponieważ w takim przypadku każdy kolejny pacjent byłby klasyfikowany jako osoba z udarem mózgu. Potrzebne jest pewne wyważenie między FP, a FN dlatego kierować się będziemy metrykami Precision - czyli ile osób zaklasyfikowanych z udarem rzeczywiście ma udar (wysoka wartość precyzji implikuje niską wartość FP) oraz Recall (wysoka wartość recall implikuje niską wartość FN), jak też ich średnią harmoniczną $F_1$.

Podsumowując, dla klasy pacjentów z udarem bardziej istotna jest metryka Recall, która informuje nas jaka jest proporcja osób zaklasyfikowanych z udarem do wszystkich osób, które zaklasyfikowane być powinny. Oczekujemy wysoką wartość, gdyż oznaczać to będzie że prawidłowo wyszukujemy osoby z udarem.

Dla pacjentów bez udaru bardziej patrzeć będziemy na metrykę Precision, która z kolei mówi nam ile zaklasyfikowanych jako osoby bez udaru rzeczywiście udaru nie mają. Wysoka wartość będzie oznaczać, że precyzyjnie wyszukujemy osoby bez udaru.

Modelowanie

powrót do spisu treści

W naszym zestawie danych występuje wiele zmiennych po wprowadzeniu polynomial features. O ile w trakcie preprocessingu danych nie usunęliśmy ich zbyt wiele i wciąż liczba zmiennych jest o wiele mniejsza od liczby obserwacji to wciąż nadmierna ich ilość prowadzić może do zjawiska overfittingu, wpływać na wyjaśnialność, czy też wprowadzać niepotrzebny szum do modelu.

Dlatego też chcemy tą liczbę zredukować poprzez wybranie naijstotniejszych zmiennych (chociaż możliwe są inne sposoby jak redukcja wymiarowości przykładowo poprzez PCA).

Wytrenujemy 3 modele:

1. Regresje logistyczną

2. SVM czyli maszyne wektorów nośnych.

3. RandomForest

Każdy z modeli zostanie wytrenowany po wcześniejszej selekcji cech z wykorzystaniem regresji logistycznej z regularyzacją LASSO. Regularyzacja ta zeruje pewne współczynniki i zostanie ona wytrenowana przy wykorzystaniu kilku parametrów regularyzujących spośród których wybrany zostanie najlepszy dla metryki ROC AUC.

Po wybraniu tych zmiennych, wytrenujemy na tym wyselekcjowanym zbiorze wszystkie modele wraz z przeszukaniem siatki hiperparametrów każdego z nich maksymalizując wynik ROC AUC co pomoże nam w doborze ostatecznego modelu.

W ostatnim kroku przejdziemy do doboru thresholdu $p$, czyli wartości prawdopodobieństwa od którego klasyfikować będziemy pacjentów z udarem. Dla regresji logistycznej oraz RFC jest to naturalne, gdyż metody te są z natury probabilistyczne a ich algorytmy zwracają wartości prawdopodobieństw klasyfikacji. Dla SVM możemy ustawić wartość hiperparametru probability=True, dzięki czemu jesteśmy w stanie pozyskać owe prawdopodobieństwa.

Całość odbywa się za scenami przy wykorzystaniu regresji logistycznej na danych wynikach SVM (algorytm Platta). Technika ta nie jest do końca wiarygodna i ma zastrzeżenia, mianowicie:

• przeprowadzana jest poprzez walidację krzyżową więc jest operacją kosztowną, zwłaszcza dla dużych zbiorów danych
• wartości prawdopodobieństw mogą nie być równe estymowanym klasom, jako że argmax wyniku nie jest równe argmax prawdopodobieństwu w każdym przypadku
• w wyniku powyższego klasa może zostać uznana jako pozytywna, mimo że prawdopodobieństwo < 0.5

W kontekście zastrzeżeń ustawiając własny threshold będziemy to w pewien sposób kontrolować tym samym ograniczając efektywność działania samego algorytmu SVM, a skupiając się na skuteczności działania techniki Platta w estymacji tych prawdopodobieństw.

Wybór techniki balansowania

powrót do spisu treści

report = pd.DataFrame()
lr_base = LogisticRegressionWithThreshold(penalty='l1', solver='saga',max_iter=10000, random_state = 123)
for X, y, name in zip([X_smote, X_smtom, X_smenn],[y_smote,y_smtom, y_smenn], ['SMOTE', 'SMOTE+TOMEK','SMOTE+ENN']):
  print('\n'+name +'\n')
  columns_to_standardize = list(set(list(X.filter(regex='age')) + list(X.filter(regex='log_glucose_no_outliers'))))

  results = fit_pred_score(X, y, X_test, y_test, lr_base, dataset_name = name, model_name = 'Base LR', scaler = True, ohe=True, columns_to_standardize = columns_to_standardize, cat_features = cat_features, visualize_test = True, visualize_train = True)
  
  report = pd.concat([report, results])
  report.reset_index(inplace=True, drop=True)

report.sort_values(by=['model_name','metric','dataset_name'])

SMOTE

SMOTE+TOMEK

SMOTE+ENN

	No stroke	Stroke	accuracy	macro avg	weighted avg	metric	dataset_name	model_name	split
2	0.840624	0.840233	0.840429	0.840429	0.840429	f1-score	SMOTE	Base LR	train
6	0.908547	0.246575	0.836892	0.577561	0.872818	f1-score	SMOTE	Base LR	test
18	0.821209	0.820244	0.820728	0.820727	0.820743	f1-score	SMOTE+ENN	Base LR	train
22	0.886301	0.265487	0.803084	0.575894	0.852794	f1-score	SMOTE+ENN	Base LR	test
10	0.842220	0.842171	0.842195	0.842195	0.842195	f1-score	SMOTE+TOMEK	Base LR	train
14	0.908364	0.258760	0.836892	0.583562	0.873302	f1-score	SMOTE+TOMEK	Base LR	test
0	0.839597	0.841265	0.840429	0.840431	0.840431	precision	SMOTE	Base LR	train
4	0.967422	0.164234	0.836892	0.565828	0.924071	precision	SMOTE	Base LR	test
16	0.847005	0.795874	0.820728	0.821440	0.822289	precision	SMOTE+ENN	Base LR	train
20	0.976604	0.166205	0.803084	0.571404	0.932863	precision	SMOTE+ENN	Base LR	test
8	0.842089	0.842302	0.842195	0.842195	0.842195	precision	SMOTE+TOMEK	Base LR	train
12	0.969417	0.171429	0.836892	0.570423	0.926346	precision	SMOTE+TOMEK	Base LR	test
1	0.841654	0.839204	0.840429	0.840429	0.840429	recall	SMOTE	Base LR	train
5	0.856426	0.494505	0.836892	0.675466	0.836892	recall	SMOTE	Base LR	test
17	0.796937	0.846154	0.820728	0.821546	0.820728	recall	SMOTE+ENN	Base LR	train
21	0.811285	0.659341	0.803084	0.735313	0.803084	recall	SMOTE+ENN	Base LR	test
9	0.842351	0.842040	0.842195	0.842195	0.842195	recall	SMOTE+TOMEK	Base LR	train
13	0.854545	0.527473	0.836892	0.691009	0.836892	recall	SMOTE+TOMEK	Base LR	test
3	3265.000000	3265.000000	0.840429	6530.000000	6530.000000	support	SMOTE	Base LR	train
7	1595.000000	91.000000	0.836892	1686.000000	1686.000000	support	SMOTE	Base LR	test
19	3265.000000	3055.000000	0.820728	6320.000000	6320.000000	support	SMOTE+ENN	Base LR	train
23	1595.000000	91.000000	0.803084	1686.000000	1686.000000	support	SMOTE+ENN	Base LR	test
11	3216.000000	3216.000000	0.842195	6432.000000	6432.000000	support	SMOTE+TOMEK	Base LR	train
15	1595.000000	91.000000	0.836892	1686.000000	1686.000000	support	SMOTE+TOMEK	Base LR	test

plot_datasets_scores = transform_report_to_plot(report[(report.model_name=='Base LR') & (report.split=='test')], x='dataset_name', cl='metric')
plotting(data = plot_datasets_scores, x='', y='value', cl='metric', type='lineplot')
print("Wykres zależności wartości wybranych metryk od techniki balansowania")

Wykres zależności wartości wybranych metryk od techniki balansowania

report_diff = prepare_diff_dataset(report, by='dataset_name')

plot_diff_scores = transform_report_to_plot(report_diff[(report_diff.model_name=='Base LR') & (report_diff.train_test_diffaccuracy.isna()==False)], x='dataset_name', cl='metric', train_test=True)
plotting(data = plot_diff_scores, x='', y='value', cl='metric', type='lineplot')
print("Wykres zależności różnic między wartościami zbioru treningowego i testowego wybranych metryk od techniki balansowania")

Wykres zależności różnic między wartościami zbioru treningowego i testowego wybranych metryk od techniki balansowania

Dla ogólnego spojrzenia spoglądamy na metrykę 'macro-avg' jako, że nie uwzględnia proporcji klas (pod kątem liczby wystąpień), a chcemy aby klasa 'Stroke' była tak samo ważna w ocenie modelu.

Tak też dla:

• metryki $F_1$ najlepszą średnią wartość uzyskał SMOTE+TOMEK
• metryki Precision najlepszą wartość dla klasy bez udaru uzyskał SMOTE+ENN
• metryki Recall najlepszą wartość dla klasy z udarem uzyskał SMOTE+ENN

W ogólności widać, że model sprawuje się całkiem dobrze bez przeprowadzenia preprocessingu oraz doboru hiperparametrów. W naszym przypadku wybór metody balansowania danych raczej nie jest zbyt istotny, ale dobrano 3 możliwości w celu pokazania różnych technik jak też pewnego wpływu na ewentualne wyniki. Tutaj różnice są naprawdę niewielkie, ale z racji na wyższe wartości, jak też niższe różnice względem zbioru treningowego i testowego wybierzemy metodą hybrydową SMOTE+ENN.

X = X_smenn.copy()
y = y_smenn.copy()
dataset_name = 'SMOTE+ENN'
columns_to_standardize = list(set(list(X.filter(regex='age')) + list(X.filter(regex='log_glucose_no_outliers')))) + ['bmi_categorical_encoded'] #bmi jest skategoryzowana, ale medianą w celu zachowania pewnej informacji w wartości dlatego ją standaryzujemy

X_train_preprocessed, X_test_preprocessed, standard_scaler = preprocess_data(X, X_test, scaler = True, ohe=True, columns_to_standardize = columns_to_standardize, cat_features = cat_features)

Selekcja zmiennych przy użyciu LASSO

powrót do spisu treści

Do selekcji zmiennych wykorzystamy metodę LASSO w tym przypadku będzie to regresja logistyczna z funkcją kary w metryce L1. Przeszukamy siatkę parametrów regulujących C (im niższa wartość tym bardziej restrykcyjna funkcja kary). Będziemy maksymalizować wynik ROC AUC, czyli maksymalizować ilość prawdziwych pozytywnych predykcji przy jednoczesnej minimalizacji negatywnych predykcji. Następnie z najlepszego estymatora wybierzemy odpowiednie zmienne przy użyciu modułu SelectFromModel.

lr_grid = {'C': [1, 0.01, 0.001, 0.0001, 0.00001]}
lr = LogisticRegressionWithThreshold(penalty='l1', solver='saga',max_iter=10000, random_state = 123)
lr_search = GridSearchCV(estimator = lr, param_grid = lr_grid, scoring = 'f1',
                          cv = 3, n_jobs = -1)
lr_search.fit(X_train_preprocessed,y)

GridSearchCV(cv=3,
             estimator=LogisticRegressionWithThreshold(max_iter=10000,
                                                       penalty='l1',
                                                       random_state=123,
                                                       solver='saga'),
             n_jobs=-1, param_grid={'C': [1, 0.01, 0.001, 0.0001, 1e-05]},
             scoring='f1')

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

lr = lr_search.best_estimator_
selection_model = SelectFromModel(lr).fit(X_train_preprocessed, y)
coefs = {}
for coef, feature in zip(selection_model.estimator.coef_[0], selection_model.estimator.feature_names_in_):
  coefs[feature] = coef

X_selected = selection_model.transform(X_train_preprocessed)
print("Wymiar początkowego zbioru danych: ", X_train_preprocessed.shape)
print("Wymiar zbioru po selekcji zmiennych: ", X_selected.shape)

coefs

Wymiar początkowego zbioru danych:  (6320, 26)
Wymiar zbioru po selekcji zmiennych:  (6320, 26)

{'heart_disease log_glucose_no_outliers': 0.8344674630857023,
 'is_rural': -0.3836178773412988,
 'work_type_encoded': 0.06631646886871749,
 'age hypertension': -0.5360502118541486,
 'is_married log_glucose_no_outliers': -1.3696022295920762,
 'is_female': -0.5486473753957081,
 'age bmi_categorical_encoded': -1.894616593787902,
 'age is_female': -0.28609178877921704,
 'age^2': -1.7626312272903806,
 'age log_glucose_no_outliers': 3.542685050407526,
 'age': 3.29947886926461,
 'is_female log_glucose_no_outliers': 0.5107077922181565,
 'is_rural log_glucose_no_outliers': 0.3124455072056057,
 'bmi_categorical_encoded': -0.27217306234852345,
 'bmi_categorical_encoded log_glucose_no_outliers': 1.213534468528651,
 'age is_rural': -0.24886293590718478,
 'log_glucose_no_outliers': -3.2697571834283563,
 'work_type_encoded log_glucose_no_outliers': 0.29329202605187144,
 'log_glucose_no_outliers^2': 2.320938764383206,
 'is_married': 4.113051514840922,
 'hypertension': -1.8745757969323773,
 'heart_disease': -1.5783786106522788,
 'age heart_disease': -0.392978351567816,
 'age is_married': -0.20652025298948645,
 'age work_type_encoded': -0.3296512727776314,
 'hypertension log_glucose_no_outliers': 1.044660005331124}

Widzimy, że żadna zmienna nie została "usuniętą" przez LASSO poprzez przypisanie zerowej wagi.

Hyperparameters tuning

powrót do spisu treści

Ostatnim etapem będzie dostrojenie hiperparametrów dla każdego z modeli. W regresji logistycznej dobierzemy wartości parametru $C$ odpowiadającego za regularyzację. Dla SVM wypróbujemy inne wartości gamma, $C$ oraz jądro liniowe. Dla RFC stworzymy przestrzeń hiperparametrów uwzględniającą m.in. głębokość drzewa. Wykorzystamy do tego GridSearchCV z pakietu sklearn. Metryką, którą będziemy optymalizować będzie metryka ROC AUC, która odzwierciedla wartości przedstawione na powyższym wykresie. Im wyższa wartość AUC tym wyższy odsetek klsayfikacji prawdziwie pozytywnych przy niskim odsetku klasyfikacji fałszywie pozytywnych. Ta metryka pozwala dobrze zestawiać ze sobą modele oraz podobnie jak poprzednie obrane - poprawnie odzwierciedla nasz cel.

svm_grid = {'C': [1, 0.01, 0.0001],
              'gamma': [1, 0.01, 0.0001],
              'kernel': ['rbf', 'linear']}
rfc_grid = {
    'bootstrap': [True],
    'max_depth': [80, 100],
    'max_features': [2, 3],
    'min_samples_leaf': [3, 5],
    'min_samples_split': [8,12],
    'n_estimators': [100, 500, 1000]
}

# Create a based model
svm = SupportVectorMachineWithThreshold()
rfc = RandomForestClassifierWithThreshold(random_state = 123)
# Instantiate the grid search model

svm_search = GridSearchCV(estimator = svm, param_grid = svm_grid, scoring = 'f1',
                          cv = 3, n_jobs = -1, verbose = 2)
rfc_search = GridSearchCV(estimator = rfc, param_grid = rfc_grid, scoring = 'f1',
                          cv = 3, n_jobs = -1, verbose = 2)

for grid in [svm_search, rfc_search]:
  grid.fit(X_selected,y)
  print(grid.best_estimator_)

Fitting 3 folds for each of 18 candidates, totalling 54 fits
SupportVectorMachineWithThreshold(C=1, gamma=1)
Fitting 3 folds for each of 48 candidates, totalling 144 fits
RandomForestClassifierWithThreshold(max_depth=80, max_features=3,
                                    min_samples_leaf=3, min_samples_split=8,
                                    random_state=123)

Otrzymaliśmy więc 3 modele z najlepszymi parametrami pod względem metryki ROC AUC z przeszukiwanych siatek parametrów. Teraz zestawimy je ze sobą.

plt.figure(0).clf()

svm = SupportVectorMachineWithThreshold(gamma=1, C=1, kernel='rbf', probability=True, random_state = 123)
rfc = RandomForestClassifierWithThreshold(max_depth=80, max_features=3,
                                    min_samples_leaf=3, min_samples_split=8,
                                    n_estimators=500, random_state=123)

#fit logistic regression model and plot ROC curve
lr.fit(X_selected, y)
svm.fit(X_selected, y)
rfc.fit(X_selected, y)

RandomForestClassifierWithThreshold(max_depth=80, max_features=3,
                                    min_samples_leaf=3, min_samples_split=8,
                                    n_estimators=500, random_state=123)

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

<Figure size 576x396 with 0 Axes>

y_pred_lr = lr.predict(X_test_preprocessed)
y_pred_svm = svm.predict(X_test_preprocessed)
y_pred_rfc = rfc.predict(X_test_preprocessed)

fpr_lr, tpr_lr, _ = metrics.roc_curve(y_test, y_pred_lr)
auc = round(metrics.roc_auc_score(y_test, y_pred_lr), 4)
plt.plot(fpr_lr,tpr_lr,label="Logistic Regression, AUC="+str(auc))

fpr_svm, tpr_svm, _ = metrics.roc_curve(y_test, y_pred_svm)
auc = round(metrics.roc_auc_score(y_test, y_pred_svm), 4)
plt.plot(fpr_svm,tpr_svm,label="Support Vector Machine, AUC="+str(auc))

fpr_rfc, tpr_rfc, _ = metrics.roc_curve(y_test, y_pred_rfc)
auc = round(metrics.roc_auc_score(y_test, y_pred_rfc), 4)
plt.plot(fpr_rfc,tpr_rfc,label="Random Forest Classifier, AUC="+str(auc))

#add legend
plt.legend()
plt.show()

Widzimy, że spośród najlepszych estymatorów dla wybranej przestrzeni parametrów najlepszy wynik osiąga regresja logistyczna. Z kolei SVM oraz RFC sprawują się niewiele lepiej jak bazowy estymator czyli predykujący z prawdopodobieństwem 50%. To tak jakbyśmy typowali czy ktoś ma udar rzutem monetą - SVM i RFC sprawują się tylko niewiele lepiej dla tego scenariusza. Ograniczymy się więc do regresji logistycznej i sprawdzimy czy wybór thresholda $p$ wpłynie na nasz wynik.

Wybór thresholdu $p$

powrót do spisu treści

for p in [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]:
    print('\nWartość p='+str(p)+' dla modelu Regresji Logistycznej\n')
    model_name = 'Logistic Regression with p='+str(p)
    results = fit_pred_score(X_selected, y, X_test_preprocessed, y_test, lr, dataset_name = dataset_name, model_name = model_name, scaler = False, ohe=False, 
                             visualize_test = True, visualize_train = True, threshold = p)
    
    report = pd.concat([report, results])
    report.reset_index(inplace=True, drop=True)

report.sort_values(by=['model_name','metric','dataset_name'])

Wartość p=0.1 dla modelu Regresji Logistycznej

Wartość p=0.2 dla modelu Regresji Logistycznej

Wartość p=0.3 dla modelu Regresji Logistycznej

Wartość p=0.4 dla modelu Regresji Logistycznej

Wartość p=0.5 dla modelu Regresji Logistycznej

Wartość p=0.6 dla modelu Regresji Logistycznej

Wartość p=0.7 dla modelu Regresji Logistycznej

Wartość p=0.8 dla modelu Regresji Logistycznej

Wartość p=0.9 dla modelu Regresji Logistycznej

	No stroke	Stroke	accuracy	macro avg	weighted avg	metric	dataset_name	model_name	split
2	0.840624	0.840233	0.840429	0.840429	0.840429	f1-score	SMOTE	Base LR	train
6	0.908547	0.246575	0.836892	0.577561	0.872818	f1-score	SMOTE	Base LR	test
18	0.821209	0.820244	0.820728	0.820727	0.820743	f1-score	SMOTE+ENN	Base LR	train
22	0.886301	0.265487	0.803084	0.575894	0.852794	f1-score	SMOTE+ENN	Base LR	test
10	0.842220	0.842171	0.842195	0.842195	0.842195	f1-score	SMOTE+TOMEK	Base LR	train
...	...	...	...	...	...	...	...	...	...
92	0.948214	0.666667	0.947212	0.807440	0.933018	precision	SMOTE+ENN	Logistic Regression with p=0.9	test
89	0.996937	0.284124	0.652373	0.640531	0.652373	recall	SMOTE+ENN	Logistic Regression with p=0.9	train
93	0.998746	0.043956	0.947212	0.521351	0.947212	recall	SMOTE+ENN	Logistic Regression with p=0.9	test
91	3265.000000	3055.000000	0.652373	6320.000000	6320.000000	support	SMOTE+ENN	Logistic Regression with p=0.9	train
95	1595.000000	91.000000	0.947212	1686.000000	1686.000000	support	SMOTE+ENN	Logistic Regression with p=0.9	test

96 rows × 9 columns

plot_datasets_scores = transform_report_to_plot(report[(report.model_name.str.contains('Logistic Regression with p')) & (report.split=='test')], x='model_name')
plotting(data = plot_datasets_scores, x='', y='value', cl='metric', type='lineplot')

0

report_diff_p = prepare_diff_dataset(report, by='model_name')

plot_diff_scores_p = transform_report_to_plot(report_diff_p[(report_diff_p.model_name.str.contains('Logistic Regression with p')) & (report_diff_p.train_test_diffaccuracy.isna()==False)], x='model_name', cl='metric', train_test=True)
plotting(data = plot_diff_scores_p, x='', y='value', cl='metric', type='lineplot')
print("Wykres zależności różnic między wartościami zbioru treningowego i testowego wybranych metryk od wartości p")
#plot_diff_scores_p

Wykres zależności różnic między wartościami zbioru treningowego i testowego wybranych metryk od wartości p

Wnioski dla:

• wartości Precision:

  • dla klasy bez udaru jest równa prawie 1 i maleje nieznacznie ze wzrostem wartości $p$. Jest to spowodowane częstym przypisaniem osób do klasy z udarem poprzez niską wartość prawdopodobieństwa, a więc ograniczenie liczby fałszywie pozytywnych dla klasy bez udaru.
  • dla klasy z udarem jest wartość niska i rosnąca wraz ze wzrostem parametru $p$. Zbiór jest niezbalansowany i przy niskim prawdopodobieństwie często będziemy klasyfikować nieprawdziwie ludzi bez udaru jako osoby z udarem.

• wartości Recall:

  • dla klasy bez udaru znacznie rośnie wraz ze wzrostem thresholdu $p$. Wynika to faktu, że im wyższy wartość $p$ tym częściej przypisujemy osoby do klasy bez udaru i minimalizujemy tym samym wartość fałszywie negatywnych dla tej klasy.
  • dla klasy z udarem również zaczyna się od wartości bliskiej 1 (podobnie jak Precyzja), ale maleje znacząco ze wzrostem wartości $p$. Jest to naturalna konsekwencja, że wraz ze wzrostem tego parametru maleje ilość przypisań do klasy z udarem (musi być wykazane wysokie prawdopodobieństwo przynależności), a więc wzrasta ilość pacjentów "pominiętych" przy przypisaniu do klasy z udarem.

Jak wspomnieliśmy na początku, chcemy maksymalizować wartość Precyzji dla klasy bez udaru (0) oraz Recall dla klasy z udarem (1). Dodatkowo wysokie wartości Precyzji dla klasy (1) oraz Recall dla klasy (0) będą atutem.

Od wartości $p=0.2$ dla LR możemy zaobserwować mocny spadek Recall dla klasy z udarem. W punkcie spadku dla klasy (0) wartość precyzji jest relatywnie wysoka, a dla klasy (1) następuje niewielki spadek względem Recall. Dodatkowo niskie wartości $p$ mają niskie różnice metryk między zbiorem treningowym a testowym. Są też zgodne z "wyczuciem dziedzinowym" tj. chcemy być bardziej uważni na możliwość wystąpienia udaru i raczej wystarczy nam już niewielkie prawdopodobieństwo wystąpienia, aby zaklasyfikować pacjenta. Co prawda wzrasta ryzyko nieprawidziwej klasyfikacji, ale w najgorszym przypadku zostaną poniesione niepotrzebne koszty (z takiej racji nie może być to też wartość zbyt niska, gdyż prawie każdego od razu będziemy klasyfikować jako osobę z udarem). Z drugiej strony zbyt wysoka wartość wymaganego prawdopodobieństwa może narazić zdrowie, a nawet życie ludzkie (widać to chociażby po drastycznym spadku Recall dla klasy 1 co oznacza, że mocno wzrosła liczba osób zaklasyfikowanych jako bez udaru, gdzie w rzeczywistości on wystąpił).

Interesują nas wartości $p \in (0.2, 0.3)$, ponieważ mają wysokie wartości metryk Precision dla klasy (0) oraz Recall dla klasy (1), które ustaliliśmy na początku. Jednocześnie są to momenty, w których wartość Precision dla (1) oraz Recall dla (0) zaczyna wzrastać wraz z $F_1$ score.

Wartość $p=0.1$ klasyfikowałaby prawie każdego jako pacjenta z udarem co jest niepożądanym efektem i do takich wniosków nie jest potrzebny żaden model. Z kolei w wartościach $p > 0.3$ widzimy mocny spadek Recall dla klasy 1, a więc tracimy pewność że wszyscy pacjencji z udarem, którzy winni być zaklasyfikowani rzeczywiście zostali przypisani do tej klasy.

Zostaje decyzja między wyższymi wartościami ustawionych apriori metryk, a pewną intuicją dziedzinową, gdyż dla $p=0.2$ również istnieje ryzyko klasyfikacji zbyt wielu "zdrowych" pacjentów. Sprawdźmy więc wartości AUC dla tych dwóch parametrów, które poinformują nas jaki scenariusz maksymalizuje TPR przy minimalnej ilości FPR.

plt.figure(0).clf()

y_pred_lr1 = lr.predict(X_test_preprocessed, threshold=0.1)
y_pred_lr2 = lr.predict(X_test_preprocessed, threshold=0.2)
y_pred_lr3 = lr.predict(X_test_preprocessed, threshold=0.3)
y_pred_lr4 = lr.predict(X_test_preprocessed, threshold=0.4)


fpr_lr1, tpr_lr1, _ = metrics.roc_curve(y_test, y_pred_lr1)
auc = round(metrics.roc_auc_score(y_test, y_pred_lr1), 4)
plt.plot(fpr_lr1,tpr_lr1,label="Logistic Regression with p=0.1, AUC="+str(auc))

fpr_lr2, tpr_lr2, _ = metrics.roc_curve(y_test, y_pred_lr2)
auc = round(metrics.roc_auc_score(y_test, y_pred_lr2), 4)
plt.plot(fpr_lr2,tpr_lr2,label="Logistic Regression with p=0.2, AUC="+str(auc))

fpr_lr3, tpr_lr3, _ = metrics.roc_curve(y_test, y_pred_lr3)
auc = round(metrics.roc_auc_score(y_test, y_pred_lr3), 4)
plt.plot(fpr_lr3,tpr_lr3,label="Logistic Regression with p=0.3, AUC="+str(auc))

fpr_lr4, tpr_lr4, _ = metrics.roc_curve(y_test, y_pred_lr4)
auc = round(metrics.roc_auc_score(y_test, y_pred_lr4), 4)
plt.plot(fpr_lr4,tpr_lr4,label="Logistic Regression with p=0.4, AUC="+str(auc))

fpr_lr, tpr_lr, _ = metrics.roc_curve(y_test, y_pred_lr)
auc = round(metrics.roc_auc_score(y_test, y_pred_lr), 4)
plt.plot(fpr_lr,tpr_lr,label="Logistic Regression with p=0.5, AUC="+str(auc))

#add legend
plt.legend()
plt.show()

Jak możemy zauważyć, potwierdziły się intuicje wypisane w poprzednim komentarzu.

Osiągnęliśmy wartość $AUC=0.7652$ dla $p=0.2$ co jest lepszym wynikiem od domyślnej wartości parametru $p=0.5$, dla której $AUC=0.7314$, więc jest pewna poprawa modelu. Widzimy również fakt, że wysokie wartości metryk wybranych apriori słusznie maksymalizują wartość AUC przy czym widocznie uwzględniają ryzyko klasyfikowania większości pacjentów jako osoby z udarem. Ewidentnie jest to dobry kierunek w przypadku niezbalansowanego zbioru jak nasz.

Dla $p=0.3$ mamy $AUC=0.7521$ co także nie jest złym wynikiem i po konsultacji ze specjalistą być może zostałby wybrany, gdyż przypisuje klasę z wyższym poziomem prawdopodobieństwa.

Widać również jak ze wzrostem parametru $p$ maleje wartość $AUC$.

Ostatecznie dobieramy wartość $p=0.2$ jako wartość naszego głównego modelu. Ostatnim etapie przygotujemy cały gotowy pipeline do predykcji udaru.

Przykładowy pipeline

powrót do spisu treści

#Save final model
filename = 'final_lr_no_smoking_status_model.sav'
joblib.dump(lr, filename)

['final_lr_no_smoking_status_model.sav']

def preprocess_final(X_new: pd.DataFrame, standard_scaler):
  """
  Function for preprocessing raw data. 
  It performs even unwanted operations because of the lack of time to prepare cleaning the code. 
  """
  df_clear = X_new
  
  ### BMI
  df_clear['bmi_no_nan'] = 28.1
  df_clear['bmi_no_outliers'] = df_clear.bmi_no_nan
  df_clear = iqr_outliers(df_clear, 'bmi_no_outliers')

  df_clear.loc[(df_clear.bmi_no_outliers < 18.5), 'bmi_categorical_encoded'] = 17.0
  df_clear.loc[(df_clear.bmi_no_outliers >= 18.5) & (df_clear.bmi_no_outliers < 24.9), 'bmi_categorical_encoded'] = 22.0
  df_clear.loc[(df_clear.bmi_no_outliers >= 24.9) & (df_clear.bmi_no_outliers < 29.9), 'bmi_categorical_encoded'] = 27.0
  df_clear.loc[(df_clear.bmi_no_outliers >= 29.9) & (df_clear.bmi_no_outliers < 34.9), 'bmi_categorical_encoded'] = 32.0
  df_clear.loc[(df_clear.bmi_no_outliers >= 34.9), 'bmi_categorical_encoded'] = 38.0

  ### GENDER
  df_clear = df_clear.loc[df_clear.gender != 'Other',:]
  df_clear.loc[df_clear.gender == 'Female', 'is_female'] = 1
  df_clear.loc[df_clear.gender == 'Male', 'is_female'] = 0
  df_clear.is_female = df_clear.is_female.astype('int')

  ### AGE
  df_clear = iqr_outliers(df=df_clear, col='age')

  ### EVER_MARRIED
  df_clear.loc[df_clear.ever_married == 'Yes', 'is_married'] = 1
  df_clear.loc[df_clear.ever_married == 'No', 'is_married'] = 0
  df_clear.is_married = df_clear.is_married.astype('int')

  ### WORK_TYPE
  le = preprocessing.LabelEncoder()
  le.fit(df_clear.work_type)
  df_clear['work_type_encoded'] = le.transform(df_clear.work_type)
  
  ### RESIDENCE_TYPE
  df_clear.loc[df_clear.Residence_type == 'Rural', 'is_rural'] = 1
  df_clear.loc[df_clear.Residence_type == 'Urban', 'is_rural'] = 0 
  df_clear.is_rural = df_clear.is_rural.astype('int')


  ### AVG_GLUCOSE_LEVEL
  df_clear['log_glucose_no_outliers'] = np.log(df_clear.avg_glucose_level)
  df_clear = iqr_outliers(df_clear, 'log_glucose_no_outliers')
  

  ### POLYNOMIAL FEATURES
  poly = preprocessing.PolynomialFeatures(2)

  #select only apriori chosen variables
  X_poly = df_clear.select_dtypes(['int32','int64','float32','float64'])
  X_poly = X_poly.loc[:,~X_poly.columns.isin(['id','smoking_status', 'bmi_no_nan','bmi','avg_glucose_level','glucose_no_outliers', 'log_bmi_no_outliers', 'bmi_no_outliers','glucose_categorical_encoded'])]
  df_poly = pd.DataFrame(poly.fit_transform(X_poly), columns = poly.get_feature_names_out(X_poly.columns))

  df_poly = df_poly[['age','age bmi_categorical_encoded','age heart_disease','age hypertension','age is_female',
                    'age is_married','age is_rural','age log_glucose_no_outliers','age work_type_encoded',
                     'age^2','bmi_categorical_encoded','bmi_categorical_encoded log_glucose_no_outliers','heart_disease','heart_disease log_glucose_no_outliers',
                     'hypertension','hypertension log_glucose_no_outliers','is_female','is_female log_glucose_no_outliers','is_married','is_married log_glucose_no_outliers',
                     'is_rural','is_rural log_glucose_no_outliers','log_glucose_no_outliers','log_glucose_no_outliers^2','work_type_encoded','work_type_encoded log_glucose_no_outliers']]
  
  cat_features = ['is_rural','is_married', 'is_female', 'hypertension', 'heart_disease', 'work_type_encoded']
  columns_to_standardize = list(set(list(df_poly.filter(regex='age')) + list(df_poly.filter(regex='log_glucose_no_outliers')))) + ['bmi_categorical_encoded'] #bmi jest skategoryzowana, ale medianą w celu zachowania pewnej informacji w wartości dlatego ją standaryzujemy

  
  ### STANDARDIZE + OHE  
  X_new_preprocessed = preprocess_data(None, df_poly, scaler = True, ohe=True, columns_to_standardize = columns_to_standardize, cat_features = cat_features, test_only=True, sc=standard_scaler)

  return X_new_preprocessed

#X_new
def predict_stroke(X : pd.DataFrame(), lr_model, standard_scaler):
  """
  Function predicts the strokes based on the given features using trained model.
  """
  #create model with the best parameters
  model = lr_model

  #preprocess data
  X_preprocessed = preprocess_final(X, standard_scaler)

  predicted_strokes = model.predict(X_preprocessed, threshold=0.3)
  
  return predicted_strokes

#Load previously trained model
loaded_model = joblib.load(filename)
loaded_model

LogisticRegressionWithThreshold(C=1, max_iter=10000, penalty='l1',
                                random_state=123, solver='saga')

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

#Predict strokes
predict_stroke(df.loc[0:5,:], loaded_model, standard_scaler)

array([0, 0, 0, 0, 0, 0])

Dodatkowe źródła

• http://cejsh.icm.edu.pl/cejsh/element/bwmeta1.element.desklight-002cb3e1-70f4-4321-a489-f4ced00e9d3b
• https://www.analyticsvidhya.com/blog/2020/10/overcoming-class-imbalance-using-smote-techniques/
• https://towardsdatascience.com/imbalanced-classification-in-python-smote-enn-method-db5db06b8d50
• https://machinelearningmastery.com/polynomial-features-transforms-for-machine-learning/
• https://towardsdatascience.com/regularization-in-machine-learning-connecting-the-dots-c6e030bfaddd
• https://en.wikipedia.org/wiki/Confusion_matrix
• https://scikit-learn.org/stable/index.html
• https://home.cs.colorado.edu/~mozer/Teaching/syllabi/6622/papers/Platt1999.pdf